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The  primary  task  addressed  ifr-tMs  report^ consists  in  representing  collocation 
results  in  terms  of  the  spherical -harmonic  expansion  of  the  geopotential.  In 
particular,  an  equilateral  grid  of  geoid  undulations  (referring  to  a  higher  order 
surface  than  an  ellipsoid)  predicted  through  the  modified  collocation  with  noise 
at  the  smoothing  level  (n'.n1)  can  be  utilized  in  a  numerical  integration 
algorithm,  eventually  producing  an  (n'.n1)  set  of  spherical -harmonic  coefficients 
consistent  with  the  collocation  results.  The  conditions  under  which  the  con¬ 
sistency  requirement  is  satisfied  are  analyzed  in  computer  simulations. 

As  the  most  important  outcome  of  these  simulations,  the  familiar  rule  of  thumb 
is  singled  out  as  an  accurate  and  clearcut  indicator  of  the  highest  degree  and 
order  spherical-harmonic  model  (N,N)  which  can  still  faithfully  represent  the 
gravity  field  as  described  by  a  discrete  set  of  data  —  here  an  equilateral  grid 
of  geoid  undulations.  This  rule  stipulates  that  N=18Q°/iL?,  where  e°  symbolizes 
the  grid  interval  in  angular  measure. ^>It  is  observed  that  up  to  the  degree 
n'=N,  the  r.m.s.  of  the  errors  in  geoTcl  unduTations  imputable  to  the  numerical 
integration  is  comfortably  small  and  increases  at  a  very  low  rate  for  increasing 
n‘,  while  beyond  the  degree  N  the  situation  worsens  drastically.  Similar  results 
are  produced  when  errors  in  gravity  anomalies  are  considered. 

Several  other  topics  are  addressed  in  this  report  as  well,  such  as  the  com¬ 
parison  between  the  single-  and  the  double-layer  point-mass  adjustment,  the 
detailed  gravity  field  resolution  in  areas  of  special  interest  via  third-phase 
adjustments,  the  possibility  of  introducing  tidal  point  masses  into  the  existing 
point-mass  adjustment,  the  utilization  of  altimeter  residuals  for  a  detection  of 
bathymetric  and  other  oceanic  anomalies,  and  the  global  representation  of  tidal 
effects  with  the  aid  of  spherical-harmonic  tidal  coefficients.  The  last  subject 
is  especially  important,  in  that  it  leads  to  the  inclusion  of  new  and  accurate 
information  into  the  adjustment  of  specific  tidal  parameters  (the  tidal  amplitude 
and  the  Greenwich  phase  angle  per  constituent)  within  the  overall  adjustment  of 
satellite  altimetry. 
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1.  INTRODUCTION 


The  initial,  or  first-phase  treatment  of  satellite  altimeter  data  has  been 
carried  out  at  AFGL  in  a  global  short-arc  adjustment  of  spherical -harmonic  (S.H.) 
potential  coefficients,  state  vector  parameters  and,  optionally,  certain  tidal 
parameters.  This  least-squares  adjustment  is  aimed  at  a  long-wavelength  resolu¬ 
tion  of  the  earth's  gravity  field  and  its  fundamental  surface,  the  geoid.  It 
has  been  the  subject  of  various  AFGL  reports  and  papers,  e.g.,  CBlaha,  1981,  19823, 
and  has  been  applied  especially  to  SEASAT  altimetry. 

For  the  reasons  recapitulated  in  CBlaha,  19813,  the  orbital  arcs  entering 
the  adjustment  have  been  limited  to  7  minutes  in  duration,  i.e.,  to  25.0°  in 
angular  length.  Most  of  the  satellite  passes  have  been  subdivided  at  the  pre¬ 
processing  level  to  satisfy  this  requirement,  except  for  natural  disruptions  in 
the  flow  of  data  (gaps  due  to  land  masses,  malfunction  or  the  altimeter,  etc.). 

Arcs  shorter  than  3°  have  been  eliminated  altogether.  Because  of  the  (14,14) 
degree  and  order  truncation  adopted  in  the  S.H.  adjustment  model,  the  shortest 
half -wavelength  of  the  geoidal  resolution  amounts  to  approximately  12.9°  (this 
is  180°/n,  with  n=14,  as  given  by  the  familiar  rule  of  thumb).  Accordingly,  the 
arc  length  of  25°  corresponds  to  the  full  wavelength  of  the  smallest  geoidal 
features  to  be  represented  by  the  adjustment. 

In  order  to  allow  the  fine  structure  of  the  earth's  gravity  field  to  be  added 
to  its  long-wavelength  features  representing  a  new  "reference  field",  second-phase 
techniques  have  been  conceived  based  on  the  results  from  the  first  phase.  One 
such  technique  has  been  designed  in  terms  of  the  point-mass  (P.M.)  parameters, 
as  described  in  several  AFGL  reports,  e.g.,  CBlaha,  1983,  19843.  Due  to  the 
banded  structure  of  the  matrix  of  normal  equations,  the  P.M.  parameters  can  be 
resolved  in  overlapping  strips,  which  can  be  combined  together  so  as  to  cover  the 
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whole  oceanic  geoid.  Optionally,  certain  tidal  parameters  can  also  be  included 
in  such  a  least-squares  adjustment;  in  this  case  the  matrix  of  normal  equations 
becomes  banded- bordered,  having  the  tidal  parameters  in  the  border. 

The  P.M.  approach  has  recently  been  complemented  by  a  specific  modification 
of  the  least-squares  collocation  with  noise,  which  provides  another  means  for  a 
second-phase  gravity  field  resolution  on  a  local,  regional,  or  global  scale.  The 
modification  resides  in  pushing  a  part  of  the  signal,  corresponding  to  the  gravity 
field  variations  beyond  the  desired  smoothing  level,  into  the  realm  of  "noise". 

A  detailed  description  of  this  technique  is  offered  in  CBlaha,  19843. 

Both  second-phase  approaches  above  are  built  on  the  altimeter  residuals 
(taken  as  minus  the  geoidal  residuals)  from  the  first  phase,  utilized  in  the  role 
of  new  observations.  Thus,  the  high-resolution  altimeter  information  enters  the 
second  phase,  while  the  state  vector  parameters,  the  adjustment  of  which  is  the 
most  time-consuming  element  in  the  treatment  of  altimeter  data,  are  no  longer 
considered  at  this  stage. 

The  above  two  approaches  are  the  topic  of  Chapter  2,  which  recapitulates  their 
basic  features  and  describes  their  ability  to  address  particular  tasks.  One  such 
task  is  concerned  with  designing  a  third-phase  approach  in  view  of  a  detailed 
gravity  field  representation  in  areas  of  special  interest.  Chapter  3  deals  with 
the  P.M.  approach  alone,  in  considering  a  double-layer  version  with  its  advantages 
and  disadvantages. 

In  addition  to  describing  the  earth's  gravity  field,  the  outcome  of  the  col¬ 
location  approach  can  also  serve  in  generating  a  set  of  S.H.  potential  coefficients 
which,  under  certain  conditions,  will  accomplish  the  same  task  equally  well.  The 
process  of  representing  collocation  results  by  a  S.H.  expansion  is  outlined  in 
Chapter  4,  while  the  conditions  under  which  such  a  representation  is  consistent 
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with  the  collocation  outcome  are  analyzed  in  Chapter  5.  This  analysis  is 
carried  out  primarily  through  computer  simulations,  in  which  extensive  use  is 
made  of  Legendre  polynomials  and  associated  Legendre  functions.  An  accurate 
and  efficient  evaluation  of  the  latter  is  the  subject  of  Appendix  1. 

The  outcome  of  Chapter  5  is  subsequently  exploited  ir  Chapter  6.  Section 
6.1  completes  the  development  of  Chapter  4,  and  thus  rounds  off  the  algorithm 
leading  from  the  collocation  results  to  a  detailed  S.H.  representation  of  the 
earth's  gravity  field.  And  Section  6.2  uses  some  of  the  principles  set  forth  in 
Chapter  5  in  developing  an  algorithm  for  a  global  representation  of  tidal  effects 
with  the  aid  of  S.H.  tidal  coefficients.  Chapters  4-6  together  with  Appendix  1 
could  be  combined  in  one  larger  unit  representing  the  main  thrust  of  the  present 
study. 

Subjects  having  a  weaker  link  to  the  principal  topics  contained  in  the  body 
of  this  report  are  relegated  to  Appendices  2  and  3.  They  deal,  respectively, 
with  the  possibility  of  introducing  tidal  point  masses  into  the  P.M.  adjustment, 
and  with  the  utilization  of  altimeter  residuals  for  a  detection  of  bathymetric, 
geomagnetic,  and  other  anomalies  in  open  oceans.  All  three  appendices,  as  well 
as  most  of  the  chapters,  are  presented  essentially  in  a  self-contained  manner, 
and  can  thus  be  read  independently  without  impairment  to  a  good  understanding 
of  the  presented  material. 

In  addition  to  .is  Final  Report,  the  reports  dealing  with  satellite 
altimetry  during  the  period  of  the  present  research  contract  have  been  CBlaha, 
1983,  19843. 


2.  EVALUATION  OF  A  LARGE-SCALE  SECOND-PHASE  GRAVITY 
FIELD  REPRESENTATION  AND  A  DETAILED  RESOLUTION 
IN  AREAS  OF  SPECIAL  INTEREST 

2.1  Point-Mass  Approach 

The  point-mass  (P.M. )  adjustment  is  built  on  the  adjusted  geoid  from 
the  first  phase  which  can  represent  a  "normal  field".  The  parameters  in  this 
method  are  the  P.M.  magnitudes  associated  with  point  masses  distributed  in  an 
equilateral  grid.  The  point  masses  can  form  a  single  layer,  in  which  case 
they  are  all  at  the  same  depth  below  the  surface  of  the  reference  ellipsoid, 
or  a  double  layer,  in  which  case  there  are  twin  point  masses  at  each  location 
of  the  grid,  separated  by  a  predetermined  vertical  distance.  The  number  of 
parameters  is  the  same  in  both  these  modes  since  a  twin  point  mass  represents 
only  one  parameter;  the  magnitudes  of  the  shallower  and  deeper  point  masses  are 
equal,  only  their  signs  differ.  The  resolution  power  depends  primarily  on  the 
grid  interval,  and  to  a  much  lesser  degree  on  their  depth  and/or  vertical 
separation  of  the  twin  point  masses.  This  is  illustrated  in  Table  1  of  the  next 
chapter.  Thus,  for  example,  if  the  point  masses  (single  or  twin)  are  distributed 
in  a  2°x2°  equilateral  grid,  the  resolution  corresponds  approximately  to  a 
(90,90)  spherical-harmonic  expansion. 

The  second-phase  P.M.  adjustment  treats  the  new  observations  (minus  the 
geoidal  residuals  from  the  first  phase),  the  P.M.  parameters  and,  optionally, 
certain  tidal  parameters  in  a  simultaneous  least-squares  process  leading  to  a 
more  detailed  resolution  of  the  earth's  gravity  field  on  the  global  oceanic 
scale.  In  view  of  the  banded  or  the  banded-bordered  structure  of  normal  equations, 
an  approximation  is  introduced  through  a  cut-off  distance  beyond  which  the  con¬ 
tribution  of  observations  to  a  given  P.M.  parameter  is  ignored.  Although  only 
information  representing  the  geoid  undulations  enters  the  adjustment  in  the 
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form  of  observed  quantities,  the  solved  for  P.M.  magnitudes  allow  the  evaluation 
of  other  functions  of  the  disturbing  potential  (with  respect  to  the  adopted 
"normal  field")  as  well,  given  the  geodetic  coordinates  on  the  earth's  surface. 

It  should  be  noted  that  due  to  the  relatively  high  degree  and  order  of  the 
adopted  reference  field,  such  as  the  (14,14)  spherical-harmonic  (S.H.)  field 
resolved  during  the  first  phase  as  compared  to  the  ellipsoidal  field,  the 
spherical  approximation  introduced  in  the  P.M.  model  is  inconsequential. 

The  final  outcome  of  the  two  adjustment  phases  with  respect  to  the  geoid 
and  the  earth's  gravity  field  consists  of  a  set  of  S.H.  potential  coefficients 
and  the  magnitudes  of  the  point  masses  introduced  at  the  predetermined  locations, 
which  enable  one  to  determine  geoid  undulations,  gravity  anomalies,  etc.,  at  any 
points  of  interest  (the  latter  usually  form  a  geographic  grid  to  be  used  for  the 
construction  of  contour  maps).  Each  such  quantity  is  an  algebraic  sum  of  two 
parts,  the  first  due  to  the  S.H.  coefficients  and  the  second  due  to  the  P.M. 
parameters.  The  numerical  values  of  such  quantities,  for  example  the  values  of 
geoid  undulations  covering  the  globe  in  a  sufficiently  dense  grid,  are  considered 
to  describe  the  earth's  gravity  field  to  within  the  desired  resolution 
characterized  by  an  equivalent  S.H.  expansion  (n',n'). 

If  the  desired  resolution  should  be  very  detailed,  practical  difficulties 
will  manifest  themselves  when  the  second-phase  adjustment  encompasses  the  whole 
globe  or  a  large  ocean  basin.  Since  the  P.M.  grid  is  increasingly  dense  in 
such  cases,  the  resolution  of  a  great  many  P.M.  parameters  becomes  impossible 
from  the  standpoints  of  computer  storage  and  run-time.  The  key  consideration 
is  that  all  the  P.M.  parameters  are  to  be  resolved  simultaneously  in  a  least- 
squares  adjustment,  not  one  P.M.  magnitude  at  a  time.  Although  the  banded  (or 
banded-bordered)  structure  of  normal  equations  has  made  a  ten-fold  increase  in 
the  number  of  parameters  possible  with  regard  to  a  general  least-squares 


algorithm,  even  this  method  will  fail  for  a  much  finer  than  2°-resolution. 
Considering  the  current  computer  capabilities,  a  1°-  or  finer  resolution  can 
be  achieved,  on  a  routine  basis,  only  in  limited  areas  of  special  interest. 

A  P.M.  adjustment  can  thus  be  conceived  in  two  steps.  One  step  constitutes 
what  has  been  called  a  second-phase  adjustment  built  on  the  global  first-phase 
adjustment.  And  the  other  step  can  be  called  a  third-phase  adjustment,  which 
amounts  to  an  additional  P.M.  densif ication  in  areas  of  special  interest.  In 
analogy  to  the  former  case,  this  adjustment  is  built  on  the  minus  geoidal  residual 
from  the  previous  step,  that  is,  from  the  second-phase  solution.  The  set  of  point 
masses  forms  now  a  denser  grid  compared  to  the  second  phase,  and  the  depth  is 
commensurably  smaller.  The  final  quantities  of  interest  are  obtained  through 
the  sum  of  three  parts,  each  corresponding  to  the  appropriate  phase.  In  this 
way,  the  transition  of  contour  lines  between  the  lower  and  the  higher  resolution 
areas  is  gradual.  As  an  example  of  SEASAT  altimetry,  the  second  phase  could  be 
designed  to  model  a  2°-geoid  over  the  global  oceans,  and  the  third  phase  could 
be  designed  to  model  a  l°-geoid  in  the  western  part  of  the  North  Atlantic.  Due 
to  SEASAT  tracks  intersecting  in  an  approximate  l°xi°  equilateral  grid,  a  resolu¬ 
tion  superior  to  a  l°-geoid  would  fail  in  principle  due  to  insufficient  data. 


2.2  Collocation  Approach 


The  basic  difference  distinguishing  the  collocation  approach  from  the  P.M. 
adjustment,  from  the  computational  point  of  view,  is  that  a  simultaneous  least- 
squares  adjustment  of  several  variables  does  not  take  place.  In  principle,  only 
one  prediction  point  is  solved  for  at  a  time.  Accordingly,  instead  of  an  inver¬ 
sion  of  a  large,  albeit  strongly  patterned  system  of  normal  equations,  inversions 
of  small  matrices  take  place,  one  inversion  per  prediction  point.  The  prediction 
points  themselves  can  be  distributed  in  an  equilateral  grid  similar  (or  identical) 
to  the  P.M.  grid  considered  previously.  In  this  case  the  resolution  power,  cor¬ 
responding  toan(n',n')  S.H.  expansion,  is  similar  in  both  approaches.  Consistent 
with  this  concept,  the  realization  of  the  (n',n')  geoid  in  the  adapted  collocation 
approach  entails  pushing  the  higher-resolution  features  (with  respect  to  the 
n',n'  S.H.  expansion)  into  the  realm  of  "noise". 

The  (new)  observations  are  again  minus  the  geoidal  residuals  from  the  first 
phase,  and  they  can  again  be  limited  to  those  located  within  a  given  spherical 
cap  from  the  desired  prediction  point.  The  predictions  of  greatest  interest  are 
those  of  geoid  undulations,  obtained  through  the  use  of  the  geoidal  covariance 
function;  predictions  of  other  quantities  such  as  gravity  anomalies  can  be  ob¬ 
tained  for  the  same  location  upon  using  the  appropriate  cross-covariance  function. 
All  the  predicted  quantities  refer  to  the  "normal  field"  represented  by  an  (n,n) 
S.H.  expansion,  where  n=14  has  usually  been  chosen. 

Since  the  number  of  observations  involved  with  one  prediction  point  can  be 
made  reasonably  small  (see,  e.g.,  weighted  averaging  of  observations  located 
close  together,  as  explained  in  Section  6.4  of  CBlaha,  19843),  one  is  faced  with 
inverting  N  small  matrices  as  opposed  to  inverting  an  N*N  system  of  normal 
equations,  where  N  is  the  number  of  prediction  points.  This  allows  for  modest 
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computer  storage  requirements  which  are  set  to  accommodate,  in  one  "loop",  a 
chosen  number  of  prediction  points  together  with  the  corresponding  vectors  and 
matrices  of  small  dimensions.  Since  the  prediction  procedure  applies  to  one 
point  at  a  time,  predictions  in  subsequent  "loops"  can  be  made  (and  stored  on  a 
magnetic  tape)  for  the  whole  oceanic  geoid  essentially  in  an  arbitrarily  dense 
grid. 

The  predictions  of  the  (n',n‘)  geoid  based  on  the  (n,n)  "normal  field"  can 
be  made  with  advantage  in  the  original  equilateral  grid  of  prediction  points, 
rather  than  in  a  finer  grid.  The  original  grid  can  later  be  densified  at  will 
upon  using  the  economical  "errorless  collocation".  Such  a  procedure  utilizes 
the  original  predictions  in  the  role  of  observations  and  computes  the  pertinent 
matrices  (MT  and  H  in  CBlaha,  1984D)  with  the  value  n',  rather  than  infinity, 
implying  that  the  expansion  of  the  covariance  function  ends  with  n'  and  that  the 
geoid  beyond  (n',n')  is  considered  to  be  "zero". 

The  densified  grid  corresponding  to  the  (n',n')  S.H.  expansion  can  be  used 
for  a  variety  of  purposes.  For  example,  it  can  serve  during  the  construction  of 
a  geoidal  contour  map  (in  this  case  it  may  be  useful  to  make  it  a  fairly  dense 
geographical  grid).  In  a  very  similar  procedure,  the  original  as  well  as  the 
densified  grid  can  be  made  to  represent  other  geophysical  quantities  besides 
geoid  undulations,  such  as  gravity  anomalies.  As  another  example,  an  equilateral 
densified  grid  can  serve  in  the  computation  of  S.H.  potential  coefficients,  via 
integral  formulas,  through  the  degree  and  order  (n',n').  A  lower  degree  and 
order  expansion,  if  desired,  can  be  obtained  simply  by  a  truncation  of  the 
expansion  just  formulated,  due  to  the  familiar  orthogonality  properties  of  spheri¬ 
cal  harmonics.  As  will  transpire  in  later  chapters,  the  original  equilateral 
grid  of  prediction  points  may  in  itself  be  sufficient  for  such  a  determination 
of  S.H.  potential  coefficients. 
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In  analogy  to  the  P.M.  approach,  a  third-phase  collocation  solution  can  be 
built  on  the  results  (minus  the  geoidal  residuals)  of  a  second-phase  solution; 
in  this  case,  n  of  the  "reference  field"  should  be  replaced  by  n'.  The  final 
quantities  of  interest  are  again  obtained  as  the  sum  of  three  parts,  etc.  How¬ 
ever,  due  to  the  practical  feature  of  solving  for  one  prediction  point  at  a  time, 
such  a  detailed  gravity  field  representation  could  be  obtained  on  a  global  scale 
already  in  the  second  phase  if  desired.  By  comparison,  in  the  P.M.  approach  the 
necessity  to  produce  a  medium-detail  resolution  in  the  second-phase  adjustment 
was  dictated  by  the  fact  that  the  most  detailed  solution  would  be  impossible  to 
obtain  on  a  global  scale,  and  that  some  global  solution  is  always  desired. 

However,  the  option  of  readjusting  certain  tidal  effects  (tidal  amplitudes 
and  phase  angles  for  selected  tidal  constituents)  is  not  available  in  the  collo¬ 
cation  approach.  Such  a  readjustment  can  be  made  in  individual  ocean  basins  in 
the  second-phase  P.M.  adjustment.  If  this  readjustment  is  desired  in  conjunc¬ 
tion  with  a  detailed  description  of  the  earth's  gravity  field,  it  could  be 
achieved  through  a  combination  of  the  P.M.  and  the  collocation  approaches.  In 
particular,  one  could  perform  a  relatively  coarse  second-phase  P.M.  adjustment, 
followed  by  a  collocation  solution  based  on  minus  the  geoidal  residuals  from 
the  former.  As  an  example,  a  global  P.M.  solution  could  be  carried  out  for  a 
4°x4°  equilateral  grid  of  point  masses  (distributed  over  the  world's  oceans), 
which  would  also  include  a  readjustment  of  the  tidal  amplitudes  and  phases  for 

the  constituents  M  ,  N  and  0  ;  in  the  case  of  SEASAT,  the  constituents 

2  2  1 

S2,  K2,  Kx  and  P1  could  not  be  properly  resolved,  see  e.g.  Chapter  5  in  CBlaha, 
19821.  The  collocation  approach  could  then  follow  with  the  prediction  points 
forming  a  1°*10  equilateral  grid,  and  result  in  a  l°-global  resolution. 
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3.  ROLE  OF  THE  POINT-MASS  CONFIGURATION 
IN  DESCRIBING  THE  GRAVITY  FIELD  DETAIL 

In  CBlaha,  19833,  the  concept  of  a  point-mass  (P.M.)  adjustment  built  on 
minus  the  geoidal  residuals  from  the  global  spherical-harmonic  (S.H.)  adjust¬ 
ment  was  extended  from  a  single-  to  a  double-layer  P.M.  configuration.  The 
latter  Involves  twin  point  masses  located  along  the  same  vertical,  i.e., 
alo-g  the  same  radius  in  spherical  approximation,  separated  by  a  pre-determined 
distance  v.  These  twin  point  masses  are  identical  in  magnitude,  but  differ  in 
sign.  The  double-layer  algorithm  has  been  tested  through  a  comparison  with  the 
same  adjustment  problem  solved  by  means  of  a  single  P.M.  layer.  The  locations 
of  point  masses  in  the  shallower  layer  have  been  made  to  coincide  with  the 
locations  in  the  single-layer  approach.  As  the  separation  v  was  being  increased, 
all  the  results  were  approaching  their  counterparts  from  the  single-layer  adjust¬ 
ment.  Finally,  when  this  separation  reached  several  earth's  radii,  the  two  sets 
of  results  became  identical.  Such  a  "depth"  of  the  second  layer  has  no  physical 
meaning,  but  represents  a  mathematical  tool  in  verifying  the  two  approaches. 

Another  test  has  been  performed  in  the  double-layer  P.M.  mode  alone  upon 
the  assumption  of  linearity  in  the  model.  This  assumption  implies  that  v  must 
be  sufficiently  small  with  respect  to  the  horizontal  separation  (s)  of  the 
point  masses  and  the  depth  (d)  of  the  shallower  layer.  In  this  context,  the 
adjustment  model  for  one  observation  and  one  twin  point  mass  as  described  on 
page  68  of  CBlaha,  19833  reads: 

0.  -  (3Cij/3R1)  v  (kM) j  ,  (3.1) 

where  0i  is  the  observed  quantity  (considered  here  as  geoid  undulation),  C_ 
is  the  coefficient  of  the  P.M.  magnitude  in  the  single-layer  mode,  R1  is  the 


radius  from  the  earth's  center  associated  with  this  layer,  v  is  the  separation 
already  defined,  and  ( kM) ^  is  the  (scaled)  magnitude  of  the  point  mass  "j". 

One  can  now  consider  n  twin  point  masses  grouped  in  the  vector  X,  in  which  case 
(3.1)  becomes 

0.  =  aX  ,  (3.2a) 

a  =  v  a'  ,  (3.2b) 

where  the  row-vector  a'  is  associated  with  the  separation  v  =  unity. 

In  the  presence  of  a  redundant  number  of  observations  one  forms  the  familiar 
observation  equations  as 

V  =  AX  +  L  ,  (3.3a) 

where  V  is  the  vector  of  residuals,  A  is  the  "design"  matrix  composed  of  the  rows 
"a"  above,  X  is  the  vector  already  defined  (containing  P.M.  magnitudes  as  param¬ 
eters  to  be  determined  from  the  adjustment),  and  L  is  the  vector  of  constant 
terms  (here  consisting  directly  of  minus  the  observed  values,  due  to  the  zero 
initial  values  of  parameters).  Similar  to  (3.2b),  one  can  write 

A  =  v  A'  ,  (3.3b) 

where  A'  consists  of  rows  a1  described  following  (3.2b).  The  least-squares 
solution  for  the  parameters  reads 

X  =  -(ATPA)'1ATPL  .  (3.4) 

In  the  case  v  =  unity  ,  one  would  similarly  have 

X'=  -(A,tPA,)“1A'tPL  .  (3.5) 

With  T.  symbolizing  variance-covariance  matrices,  due  to  (3.3b)  it  follows  that 


Equation  (3,6b)  implies  that 


°x  =  (1/v)ox-  *  (3.7) 

j  j 

These  results  have  been  verified  upon  comparing  the  outcome  of  a  given  adjustment 
where  s  =  222  km,  d  =  175  km,  and  where  v  has  varied  as  v  =  500  m,  v  =  100  m, 
and  v  =  1  m. 

The  adjusted  values  of  observed  quantities  are  confirmed  to  be 

La  =  AX  =  v  A'  (1/v)  X'  =  A' X*  =  l,a  .  (3.8) 

This  implies,  under  the  assumption  of  small  v,  that  the  adjusted  observations  are 
invariable  of  v.  Similar  outcome  is  reached  for  adjusted  linear  functions  of 
parameters : 

La  =  AX  =  v  A'  (1/v)  X'  =  A'X'  =  L,a.  (3.9) 

Equations  (3.8)  and  (3.9)  have  also  been  confirmed  through  the  adjustment  mentioned 
at  the  end  of  the  last  paragraph.  The  quantities  in  (3.8)  were  geoid  undulations 
while  the  quantities  in  (3.9)  were  gravity  anomalies  and  deflections  of  the 
vertical.  In  fact,  selected  variance-covariances  of  all  these  quantities  have 
been  verified  as  well.  In  summary,  the  final  results  (geoid  undulations,  gravity 
anomalies,  etc.)  as  well  as  their  variance-covariances  are  invariant  of  the 
(small)  numerical  value  of  v,  unlike  in  the  case  of  the  P.M.  parameters  them¬ 
selves,  which  vary  in  proportion  to  1/v  and  whose  variance-covariances  vary  in 
2 

proportion  to  1/v  . 


After  being  verified  in  the  manner  described  above,  the  double-layer  P.M. 
algorithm  has  been  utilized  to  conduct  a  series  of  tests.  The  purpose  of  these 
tests  has  been  to  determine  what  kind  of  influence  can  be  exercised  by  the  value 
of  v  (no  longer  small)  on  the  fit  of  the  adjusted  P.M.  geoid  to  the  altimeter 
residuals  from  the  S.H.  adjustment.  The  measure  of  the  fit  is  provided  by  the 
root  mean  square  (r.m.s.)  of  the  resulting,  or  second-phase,  residuals.  The 
oceanic  area  selected  for  these  tests  is  near  Antarctica,  delimited  by  the 
parallels  -68°  and  -44°,  and  by  the  meridians  150°  and  280°.  This  area  produced 
22,819  residuals.  The  point  masses  were  distributed  in  a  2°x2°  equilateral  grid 
(approximately  222  km  x  222  km),  extending  beyond  the  residual  area,  so  that 
s  =  222  km.  The  number  of  (twin)  point  masses  was  677.  The  radius  of  the 
spherical  cap  beyond  which  the  observations  are  ignored  in  conjunction  with  a 
given  point  mass  was  stipulated  in  accordance  with  previous  conventions  as  1.5  s. 
The  depth  of  the  shallower  layer  was  chosen  according  to  the  d/s  ratios: 

0.1/1,  0.2/1,  0.4/1,  0.6/1,  0.8/1,  1.2/1,  and  1.6/1.  The  vertical  separation 
between  the  two  layers  was  chosen  as  v  -*-0  (in  practice  100  m) ,  0.4  d,  0.6  d, 

O. 9  d,  1.5  d,  and  v  -*■  «>  .  The  last  case  amounts  to  performing  a  single-layer 

P. M.  adjustment. 

The  above  test  cases  are  grouped  in  Table  1  featuring  their  r.m.s.  residual. 
Underneath  each  r.m.s.  value  is  listed  in  parentheses  the  average  magnitude  of 
the  residuals  (i.e.,  the  average  absolute  value).  The  average  residual  in  each 
category  approaches  zero  and,  therefore,  is  not  listed.  The  largest  values 
would  be  +0.04  m  in  the  first  row  with  the  exception  of  0.00  m  for  v  -►  °°  .  In 
all  the  other  cases  the  size  of  the  average  residual  is  2  cm  or  less.  One 
notices  that  with  the  exception  of  most  cases  in  the  categories  d/s  =  0.1/1 
and  d/s  =  0.2/1  ,  the  r.m.s.  residual  ranges  within  narrow  limits,  from  0.91  m 
to  1.00  m.  This  accords  very  well  with  the  theoretical  sigma  of  approximately 
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v=0.4d 


v=0.6d 


v=0. 9d 


v=l. 5d 


y-wo 


v-*0 


1.22 

(0.89) 

1.20 

(0.87) 

1.20 

(0.87) 

1.20 

(0.87) 

1.14 

(0.82) 

0.95 

(0.67 

1.10 

(0.78) 

1.05 

(0.74) 

1.03 

(0.75) 

1.01 

(0.75) 

0.98 

(0.69) 

0.95 

(0.67 

0.94 

(0.66) 

0.93 

(0.65) 

0.92 

(0.65) 

0.92 

(0.64) 

0.92 

(0.64) 

0.97 

(0.69 

0.91 

(0.64) 

0.91 

(0.64) 

0.92 

(0.64) 

0.92 

(0.65) 

0.93 

(0.66) 

0.  98 
(0.7C 

0.92 

(0.64) 

0.93 

(0.66) 

0.94 

(0.67) 

0.95 

(0.67) 

0.96 

(0.68) 

0.95 

(0.7C 

0.96 

(0.68) 

0.97 

(0.69) 

0.98 

(0.70) 

0.98 

(0.70) 

0.98 

(0.70) 

0.99 

(0.71 

0.98 

(0.70) 

0.99 

(0.70) 

0.99 

(0.70) 

0.99 

(0.70; 

0.99 

(0.71) 

Table  1 

f  the  test  cases  with  varying  d,  v;  listed  are  r.m.s.  of  the 
id,  in  parentheses,  their  average  magnitude  (both  in  meters) 
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one  meter,  obtained  from  the  covariance  function  for  the  corresponding  (90,90) 
S.H.  expansion.  Likewise,  the  average  magnitude  ranges  within  narrow  limits 
from  0.64  m  to  0.71  m. 

A  global  P.M.  adjustment  has  recently  been  performed  at  AFGL  in  the  single 
P.M.  mode  (v  -*■«>),  with  d/s  =  0.8/1.  It  follows  from  the  above  experience  that 
if  the  same  adjustment  were  performed  in  the  configuration  d/s  =  0.6/1  in  con¬ 
junction  with  a  small  or  moderate  value  of  v  (between  zero  and  50  km  or  even 
100  km),  the  r.m.s.  residual  would  likely  improve  by  some  7%.  This  improvement 
is  by  no  means  substantial  and  might  not  even  justify  a  global  adjustment  in  the 
double-layer  P.M.  mode  which  is  computationally  more  demanding  than  the  single¬ 
layer  mode.  Indeed,  the  main  outcome  of  the  above  tests  points  to  rather 
insignificant  variations  in  the  geoidal  fit  as  a  function  of  the  depth  of  the 
shallower  P.M.  layer  (provided  it  is  within  reasonable  limits)  and  as  a  function 
of  the  separation  between  the  two  layers.  This  indicates  that  the  power  of  the 
resolution  is  linked  almost  entirely  to  the  horizontal  distribution  of  the  point 
masses. 
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4.  REPRESENTATION  OF  COLLOCATION  RESULTS 
BY  A  SPHERICAL-HARMONIC  EXPANSION 


This  chapter  is  concerned  with  designing  a  strategy  which  will  lead  to 
compatible  gravity  field  representations  using  two  methods.  The  first  method 
is  the  modified  collocation  with  noise  as  developed  in  CBlaha,  1984D,  and  the 
second  method  comprises  a  spherical-harmonic  (S.H.)  expansion.  The  above  term 
"modified"  does  not  mean  that  the  general  concept  of  the  least-squares  collo¬ 
cation  with  noise  is  violated  in  any  sense.  Rather,  this  concept  is  applied  in 
a  modified  situation,  where  the  earth's  gravity  field  is  considered  to  be 
idealized  (smoothed  out)  in  that  it  is  exactly  expressible  by  a  S.H.  expansion 
complete  through  a  desired  degree  and  order  (n',n'). 

The  important  characteristic  of  the  second  method  is  that  the  S.H.  potential 
coefficients  are  to  be  computed  from  the  results  of  the  first  method,  which  they 
should  ideally  reproduce.  Here  these  coefficients  will  be  derived  from  the 
collocation  predictions  of  geoid  undulations.  A  similar  procedure  could,  of 
course,  be  undertaken  with  respect  to  gravity  anomalies.  But  since  the  present 
treatment  is  based  on  altimetric  observations,  the  collocation  predictions  of 
geoid  undulations  are  more  directly  derived  quantities  than  gravity  anomalies 
or  other  functions  of  the  geopotential,  and  as  such  will  be  used  in  the  analysis. 

The  general  formulas  for  the  least-squares  collocation  with  noise  read 
P  =  Mt(H  +  I)_1F  ,  (4.1a) 

O  =  H  -  M7(H  +  Z)-1M  ,  (4.1b) 

p  p 
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where 

P  =  the  vector  of  predicted  quantities,  here  geoid  undulations  at 
the  chosen  prediction  points, 

C~  =  the  error  measure  of  these  quantities,  similar  in  character  to 
a  variance-covariance  matrix, 

F  =  the  vector  of  actual  observations,  here  geoid  undulations  at  the 

given  observation  points;  it  is  composed  of  the  vector  F  (errorless 
geoid  undulations)  and  the  vector  e  (observational  noise),  namely 

F  =  F  +  e  ,  (4.1c) 

E  =  the  variance-covariance  matrix  of  the  observational  noise,  and 

H,  Mt,  Hp  =  matrices  formed  in  a  standard  way,  here  using  the  covariance 
function  for  geoid  undulations,  where  H  pertains  to  observation 
points,  Mt  to  both  prediction  and  observation  points,  and  Hp  to 
prediction  points. 

The  general  covariance  function  for  geoid  undulations  is  given  by 

oo 

D( ip)  =  E  d  P  (cos4>)  ,  (4.2) 

k=n+l  k  K 

where  n  is  the  degree  and  order  of  the  reference  field,  here  a  (14,14)  field, 

hence  n=14;  it  is  thus  apparent  that  the  term  "geoid  undulations"  as  used  above 

should  be  understood  with  respect  to  such  a  reference  field.  Further,  dk  is 

the  k-th  degree  variance  and  P  (cosip)  is  the  Legendre  polynomial  in  the  argument 

k 

cosp,  where  ip  is  the  spherical  distance  between  the  points  concerned.  Equation 
(4.2)  corresponds  to  the  actual  gravity  field,  which  could  be  thought  of  as 
described  by  an  (  oo^oo  )  S.H.  expansion. 

Suppose  now  that  an  idealized  gravity  field  is  exactly  expressible  through 
an  ( n 1 ,n * )  S.H.  expansion;  i.e.,  the  geoid  undulations  corresponding  to  the  S.H. 


coefficients  beyond  the  degree  and  order  n'  are  identically  zero.  In  such  a 
case,  the  geoidal  covariance  function  would  become 

n1 

D’M  =  E  d  P  (cost)  .  (4.3) 

k=n+l 

The  formulas  (4.1a,b)  would  remain  unchanged  in  this  situation,  except  that  the 

T 

matrices  H,  M  and  Hp  would  be  formed  with  the  aid  of  the  covariance  function 
given  by  (4.3)  instead  of  (4.2).  Furthermore,  the  vector  F  (measured  via  r) 
would  become  F' ,  composed  of  errorless  geoid  undulations  associated  with  the 
S.H.  expansion  (n',n'),  currently  under  consideration,  instead  of  the  expansion 
(°°.00) . 

Next,  make  several  temporary  assumptions  based  on  the  notion  just  introduced 

1)  the  measured  surface  is  the  idealized  geoid  corresponding  to  the 
above  gravity  field  (n',n'); 

2)  the  observations  have  an  ideal  global  configuration,  i.e.,  they  are 
sufficiently  dense  and  uniform; 

3)  the  observations  have  an  ideal  quality,  i.e.,  e  =  0,  Z  =  0; 

4)  the  predictions  are  made  using  all  of  the  observations  in  conjunction 
with  any  one  prediction  point; 

5)  the  predictions  have  an  ideal  global  configuration,  i.e.,  are  made 
in  an  arbitrarily  dense  equilateral  grid. 

Clearly,  such  assumptions  are  unrealistic  both  with  regard  to  the  observations 
(items  2  and  3)  and  to  the  computation  of  the  predictions  (items  4  and  5). 
Although  the  basic  assumption  (item  1)  is  unrealistic  as  well,  it  expresses  an 
approximation  of  the  actual  gravity  field  by  the  (n',n')  field,  which  is  an 
important  goal  common  to  both  methods.  The  other  assumptions  (items  2-5) 
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AC 

nm 

=  Cl/(4TfR)T  //  N  P  (si n<p ) 

cos  mX 

AS 

nm 

ran 

0 

sin  mX 

reflect  an  initial  stage,  where  the  two  methods  lead  to  identical  results. 

These  assumptions  will  eventually  be  relaxed,  but  one  should  keep  in  mind  that 
the  gravity  field  representation  should  be  consistent  in  the  two  methods. 

The  first  method  presents  the  idealized  geoid  already  as  the  result  of  item 
5  above,  in  that  a  geoidal  contour  map  can  be  constructed  from  the  (dense) 
global  predictions.  On  the  other  hand,  due  to  the  same  item  the  S.H.  potential 
coefficients  in  the  second  method  can  be  computed  by  the  integral  formula 


(4.4) 


where,  on  the  right-hand  side, 

R  =  the  earth's  mean  radius, 

N  =  the  geoid  undulation  (referring  to  the  14,14  reference  field), 

Pnm(sin<i))  =  the  normalized  Legendre  functions  in  the  argument  sin4>  , 

<f> , X  =  the  geocentric  latitude  and  longitude,  respectively,  of  the 
point  associated  with  N  ,  and 

do  =  the  (spherical)  surface  element  associated  with  N. 

Due  to  the  relatively  high-degree  and  order  reference  field,  the  spherical 
approximation  is  inconsequential.  For  the  same  reason,  one  has  AC  =  C  , 

nm  ran 

while  AS  S  ,  the  normalized  S.H.  coefficients, 
nm  nm 

Upon  using  the  S.H.  potential  coefficients  as  computed  above,  the  second 
method  yields  the  geoid  undulations  according  to  the  formula  which,  in  the 
present  case,  reads 

n 1  k  _ 


N  =  R  A  £  (AC,  ,  cos  iX  +  AS  sin  £X)  P  (si ncp )  , 
k=n+l  2=0 


(4.5) 
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where  4>,A  refer  to  the  point  associated  with  N.  If  such  a  point  coincides 
with  a  collocation  prediction  point,  the  predicted  value  should  be  exactly 
reproduced  by  (4.5).  Thus,  under  the  above  assumptions  the  two  methods  are 
equivalent. 

Computer  simulations  in  the  next  chapter  will  resolve  the  practical  question 
regarding  the  density  of  prediction  points  needed  for  a  satisfactory  S.H. 
representation  of  the  collocation  results.  In  this  task,  the  assumptions  (1)  - 
(4)  are  left  intact,  but  the  predictions  are  generated  in  a  finite  equilateral 
grid.  For  example,  if  n'=90,  such  a  grid  could  have  no  more  than  2.23°  on  the 
side,  due  to  a  relationship  between  the  number  of  measurements  and  the  number 
of  unknowns.  But  from  the  standpoint  of  a  good  S.H.  representation,  this  value 
should  be  smaller.  However,  an  answer  to  the  question  whether  a  20*2°  or,  per¬ 
haps,  a  l°xl°  equilateral  grid  is  necessary  in  this  task  will  be  obtained  from 
computer  simulations  which  will  generate  the  required  observations  (items  2  and 
3)  describing  the  idealized  geoid  (item  1).  In  particjiar,  one  will  be  able  to 
find  the  threshhold  size  beyond  which  the  familiar  orthogonality  relations  for 
spherical  harmonics  are  disturbed  in  a  way  significantly  affecting  the  results. 

Such  a  threshhold  can  then  serve  in  computing  the  S.H.  potential  coef¬ 
ficients  using  the  actual  altimetric  observations.  The  remaining  assumptions 
(items  1-4)  are  dealt  with  at  the  level  of  the  collocation  approach  alone.  By 
virtue  of  the  relatively  high-degree  and  order  reference  surface,  the  covariance 
function  decreases  rapidly  with  increasing  distance,  which  allows  the  use  of 
only  those  observations  which  are  located  in  a  relatively  small  spherical  cap 
centered  at  a  given  prediction  point  --  not  all  the  available  observations. 

This  consideration  removes  the  need  for  the  assumption  of  item  4  in  practical 


terms. 


The  most  important  assumption  is  that  of  item  1.  It  is  treated  in  the 
modified  collocation  approach  by  splitting  F  into  F' ,  corresponding  to  the 
(n',n‘)  field,  and  F",  representing  the  effects  beyond  this  field.  Since  the 
actual  observations  sense  the  effects  of  the  entire  field  --  not  only  those  of 
the  hypothetical  (n',n')  field  —  the  part  F"  acts  as  if  it  were  augmenting 
the  level  of  the  noise,  in  the  sense  of  adding  a  part  of  the  signal  to  the  noise. 
In  considering  (4.1c),  this  concept  is  represented  by  the  new  formula 

F  =  F'  +  (F”  +  e)  , 

where  F'  corresponds  to  the  idealized  geoid  as  just  stated,  and  (F"  +  e)  is  the 
modified  noise.  The  resulting  modified  collocation  formulation  is  similar  to 
(4.1a,b),  except  that  the  covariance  function  involved  in  the  formation  of  the 
matrices  MT  and  Hp,  but  not  H,  is  utilized  in  conjunction  with  n1  (equation  4.3) 
instead  of  °°  (equation  4.2).  Of  the  remaining  assumptions,  item  (2)  is  fulfilled 
reasonably  well  over  the  oceanic  areas,  and  item  (3)  is  inconsequential,  con¬ 
sidering  the  precision  of  satellite  altimetry. 


5.  EVALUATION  OF  SPHERICAL-HARMONIC  POTENTIAL  COEFFICIENTS 
VIA  INTEGRAL  FORMULAS 

5. 1  Design  of  the  Basic  Procedure 

The  main  topic  to  be  addressed  concerns  the  evaluation  of  spherical- 
harmonic  potential  coefficients  via  integral  formulas,  with  geoid  undulations 
as  the  known  quantities.  The  emphasis  on  geoid  undulations  stems  from  their 
availability  over  most  of  the  earth's  surface,  directly  attributable  to  an 
increased  exploitation  of  satellite  altimetry.  The  integral  formulas  will  be 
applied  to  discrete  data  sets  and  their  acceptability  will  be  assessed  through 
computer  simulations. 

The  quantities  of  key  importance  in  this  task  consist  of  "errorless"  geoid 
undulations  (N)  generated  through  the  degree  and  order  n ' ;  the  reference  ellipsoid 
is  assumed  to  have  the  same  potential  as  the  geoid  and  the  same  mass  as  the  earth, 
and  is  centered  at  the  earth's  center  of  mass.  In  the  spherical  approximation, 
the  formula  used  in  generating  geoid  undulations  reads 

n '  n  _  _  _ 

N  =  R  E  Z  (AC  cos  mX  +  AS  sin  mX)P  (sincp)  ,  (5.1) 

_  „  run  nm  ran 

n=2  m=0 

where  R  is  the  earth's  mean  radius  (6371  km),  <p  and  X  are  the  geocentric  latitude 
and  longitude,  respectively,  AC  and  AS  5  AS  are  the  normalized  spherical- 

nm  ran  ran 

harmonic  (S.H.)  potential  coefficients  referring  to  the  normal  (ellipsoidal) 

field,  and  P  (sin<j>)  are  the  normalized  associated  Legendre  functions  in  the 
ran 

argument  sin4>.  Due  to  the  potential  and  mass  assumptions,  the  mean  value  of 
geoid  undulations  is  zero. 


The  integral  formula  giving  the  S.H.  coefficients  based  on  the  knowledge 


of  N  is  written  as 


AC 

nm 

=  Cl/(4mR)3 N  P  (sin*) 

cos  mX 

AS 

nm 

nm 

0 

sin  mx 

(5.2) 


where  a  denotes  the  surface  of  the  unit  sphere  and  do  is  the  surface  element. 
This  formula  follows  from  (5.1),  where  the  geoid  undulations  are  regarded  as 
functions  on  a  sphere  and  are  expressed  in  terms  of  surface  spherical  harmonics 
(see  Sections  1-13  and  1-14  in  CHeiskanen  and  Moritz,  19673). 

Suppose  now  that  the  elements  do  are  replaced  by  finite  surface  areas  doi> 
each  associated  with  a  centrally  located  geoid  undulation  fT.  If  these  areas 
are  sufficiently  small,  the  integration  operator  in  (5.2)  can  be  replaced  by  the 
summation  operator.  If,  moreover,  all  of  doi  have  an  equal  area,  the  formula  is 
further  simplified  to  read 


AC 

nm 

p  _ 

cos  mX. 

1 

AS 

nm 

=  ( l/R) ( 1/P )  1  N.  P  (sin*.) 

i-l  1 

sin  mX. 

1 

where  P  is  the  total  number  of  geoid  undulations  fT,  whose  locations  (points  "i") 
are  described  by  the  coordinates  (<J>  ,X  ).  In  fact,  the  equal  area  criterion 
leads  to  all  of  the  Ni  in  (5.3)  having  an  equal  weight.  Thus,  in  this  "finite" 
case,  the  integral  formula  (5.2)  becomes  the  numerical  integration  formula  (5.3). 

For  the  sake  of  simplification  in  practical  computations,  the  points  "i", 
i  =  1,  2,  ...,  P,  are  usually  desired  in  a  regular  grid.  Due  to  the  difficulty 
in  computing  the  associated  Legendre  functions,  it  is  expedient  to  design  the 
grid  in  4-  and  X-  coordinates  to  that  these  functions  are  computed  only  once 


for  a  number  of  grid  points  (along  the  same  parallel).  The  grid  intervals  in 
4>  and  X  are  taken  as 


A<p  ,  A4>/cos4>i  ,  (5.4) 

respectively.  The  surface  areas  doi  are  then  considered  to  form  geographic  com¬ 
partments,  called  blocks,  whose  dimensions  A$  and  AX  are  again  given  by  (5.4), 
and  which  contain  the  points  "i"  in  their  geographic  centers.  By  virtue  of  (5.4), 
the  equal  area  requirement  for  the  blocks  is  satisfied  and  the  grid  formed  by 
the  points  "1"  is  called  equilateral.  Equation  (5.3)  can  now  be  presented  in  a 
computationally  advantageous  form: 


AC 


ran 


AS 


ran 


K 

( 1/R) (1/P)  Z  CP  (sin*)  E  N 


k=l 


ran'  k 


£=1 


k£ 


cos  mX 


k£ 


sin  mX 


k£ 


(5.5a) 


P  *  Z  P,  , 

k=l  K 


(5.5b) 


where  K  is  the  number  of  parallels  containing  grid  points  and  pk  is  the  number 
of  grid  points  on  the  parallel  "k". 

With  the  S.H.  coefficients  computed  from  (5.5a,b),  one  should  ideally 
recover  the  errorless  values  of  N  by  a  new  application  of  (5.1),  whether  at 
the  grid  points  where  N  (NkJl  in  the  notations  of  5.5a)  were  originally 
generated,  or  at  other  points.  If  the  recovery  is  judged  acceptable,  so  is 
the  density  of  the  grid.  Otherwise,  the  grid  is  too  coarse  with  respect  to 
the  degree  and  order  of  the  desired  S.H.  coefficients,  in  the  sense  that  the 
familiar  orthogonal  relations  are  no  longer  fulfilled  in  the  finite  case.  These 
relations,  demonstrated  in  Section  1-14  of  CHeiskanen  and  Moritz,  1967],  read 


JJ  R  R  do  =  0,  //  S  S  dc  =  0  ,  //  R  S  do  =  0 

nm  sr  nm  nr  "  nm  sr 


if  s^n  or  r^m  or  both 


in  any  case 


//  R2  do  =  4tt  , 
nm 

0 


//  S2  do  =  4tt  ; 
a  nm 


R 

nm 

cos  mX 

S 

nm 

=  P  (sin<|>) 

nm 

sin  mX 

Possible  difficulties  arising  from  the  subdivision  of  a  sphere  into  P 
compartments  of  finite  dimensions  can  be  assessed  by  exposing  the  connection 
between  the  equivalent  formulas  (5.1)  and  (5.2).  First,  (5.1)  is  rewritten 
with  the  notations  from  (5.6c)  as 


n  11  —  _ 

N  =  R  I  I  (AC  R  +  AS  S  ) 

-  _  nm  nm  nm  nm 

n=2  m=0 


From  here  it  follows  that 


R  ,  R 

sr  n '  n  _  _  sr 

//  N  do  =  R  I  I  (AC  JJ  R 

a  r  n=2  m=0  a  ■? 


do  +  AS  //  S 


nm  nm 


Upon  considering  (5.6a ,b),  one  readily  obtains 


Thus,  equation  (5.2)  has  been  rederived,  contingent  precisely  on  the  validity 
of  all  the  orthogonal  relations  in  (5.6a,b).  If  these  relations  are  no  longer 
valid  due  to  the  finite  dimensions  of  the  blocks  implicated  in  (5.3)  or  (5.5a,b), 
the  S.H.  coefficients  computed  by  the  latter  become  distorted  and,  consequently, 
the  errorless  values  as  obtained  in  (5.1)  with  the  original  coefficients  cannot 
be  restored  by  applying  this  formula  in  conjunction  with  the  new  coefficients. 

The  larger  the  blocks,  the  greater  the  distortions  in  the  coefficients  and  the 
greater  the  distortions  in  geoid  undulations  (differences  between  the  errorless 
and  the  recomputed  values  of  N). 


5.2  Initial  Test  Relations 


The  main  thrust  of  the  final  analysis  will  be  directed  toward  the 
numerical  assessment  of  the  distortions  in  geoid  undulations.  In  particular, 
the  root  mean  square  (r.m.s.)  value  of  these  distortions  will  be  evaluated 
against  the  background  of  a  pre-determined  standard  of  comparison  in  order  to 
judge  the  grid's  acceptability.  However,  one  may  wish  to  perform  preliminary 
tests  involving  solely  the  errorless  values  The  reasons  may  vary,  such 
as  obtaining  an  early  indication  of  an  improper  grid  design  (including  its 
density)  or  verifying  the  data-generating  program  and  its  focal  point,  the 
evaluation  of  the  associated  Legendre  functions.  Two  kinds  of  tests  are 
designed  at  this  stage.  The  first,  very  simple  test  concerns  the  average  value 
of  N.  in  conjunction  with  any  degree  of  truncation  (n‘).  And  the  second  test 
concerns  the  r.m.s.  values  of  in  two  different  truncations,  together  with 
the  r.m.s.  difference  between  these  two  sets.  The  results  of  the  second  test 
hinge  again  on  the  fulfillment  of  the  orthogonal  relations  (5.6a,b)  in  the 
finite  case. 

The  first  test  mentioned  above  is  based  on  the  familiar  relations 

{S  d0  *  SS  da  .  0  .  (5.10) 

Accordingly,  (5.7)  readily  yields 

Nq  =  (1/4tt )//  N  do  =  0  ,  (5.11) 

o 

confirming  the  earlier  statement  about  the  mean  value  of  geoid  undulations.  In 
the  present  finite  context,  (5.11)  corresponds  to 

p 

Ave(N  )  e  (1/P)  E  N.  =  0  ,  (5.12) 

i=l 

which  has  been  confirmed  with  all  the  simulations  performed  in  this  study. 


mean  value  of  CN(n")  -  N(n' )3  ,  where  the  symbols  (n‘)  and  (n“)  indicate  the 
degrees  of  truncation.  Thus,  N(n')  is  given  by  (5.7)  as  it  stands,  while 


N(n")-N(n')  can  be  written  as 


n"  n  _  _  _  _ 

N(n")  -  N(n')  =  R  Z  Z  (AC  R  +  AS  S  )  . 


n=n'+l  m=0 


nm  nm  nm  nm 


(5.13) 


Since  the  form  (5.13)  differs  from  that  of  (5.7)  merely  by  the  degrees  implied 
by  the  first  summation,  derivations  involving  any  such  expressions  follow  the 
same  pattern.  We  now  form  the  mean  of  (5.13)  over  the  unit  sphere,  namely 


n  n 

( 1/ 4tt)  Jj  CN(n")-N(n'  )32do  =  (1/4tt)R2C  Z  Z  (AC2  JTr2  do+AS2  Xfs2  do) 


n=n'+l  m=0  o 


run  nm  nm  run 


+  ...  JJ  (other  products)doD 


The  application  of  (5.6a,b)  yields 


(1/4tt)  JJ  CN(n")-N(n'  )32do  =  R2  Z  Z  (AC2  +  AS2  )  . 
0  n=n'+l  m=0  ran  nm 


(5.14) 


Since  it  can  be  shown  (see  page  91  of  CBlaha,  19843)  that 


R2  Z  (AC2  +  AS 2  )  =  d  , 
„  nm  nm  n 
m=0 


(5.15) 


where  dn  is  the  n-th  degree  variance  for  geoid  undulations,  one  can  finally  write 


(l/47T)J/CN(n")-N(n’)32do  h  M|CN(n,,)-N(n')]2|  =  Z  d  , 


(5.16) 


n=n* + 1 


where  M  designates  the  mean  value  operator.  In  analogy  to  (5.16),  we  also  have 


n'  n" 

M{[N(n' )rf  =  I  d  ,  M|CN(nM)Dzf  =  I  dn  ,  etc. 

n=2  n  n=2 

In  assuming  that  n'  <  n"  <  n  ,  a  further  application  of  (5.16)  yields 

(5.17a) 
(5.17b) 


M{CN(n,,)-N(n'):2}  +  M{CN(n)-N(n")D2i  =  M{cN(n)-N(n'  )D2  } 


n=n'+l 


n 

I 

n=n"+l 


n 

l 

n=n' +1 


the  three  terms  in  (5.17b)  are  the  equivalent  expressions  of  the  terms  directly 
above  them. 

In  the  finite  context,  the  usual  replacement  of  operators  means  that  (5.17a) 
is  replaced  by  the  mean  square  (m.s.)  relationship 

m.s.CN(n")-N(n'  )]  +  m.  s.  CN(rf)-N( n" )]  =  m. s.  CN("n)-N(n'  )D  .  (5.18a) 


In  the  sequel,  the  r.m.s.  values  are  symbolized  by  o  ,  so  that  the  m.s.  values 

-  2 

can  be  written  as  a  .  Furthermore,  the  identification  of  various  degrees  and 
their  relationships  is  made  via  lower  indices  as  is  illustrated  in  the  following 
transcription  of  (5.18a): 


-2 

n' +1 ,n" 


+  0 


n"  +1  ,TT 


'2 

°n'+l,7T 


(5.18b) 


The  usefulness  of  the  indices  transpires  especially  upon  comparing  (5.18b)  with 
(5.17b).  As  has  been  the  case  several  times  before,  the  validity  of  the  expres¬ 
sions  such  as  (5.18b)  hinges  on  the  fulfillment  of  (5.6a,b)  in  the  finite 


context. 


The  above  discussion  can  be  illustrated  by  the  following  example,  wherein 
two  sets  of  geoid  undulations  are  generated  in  a  2°x2c  equilateral  grid  contain¬ 
ing  10,328  points.  These  sets  are  identified  by  the  truncation  degrees  25  and 
70,  implying  that  n'  =2,  n"  =25,  and  n  = '/0.  The  r.m.s.  values  of  Ni  are  computed 
as 

a  =  30.2528m  ,  a  =  30.3000m  . 

4  /  4  v  c.  $  /  U 


The  differences  in  geoid  undulations  are  evaluated  at  all  the  grid  points, 
resulting  in  the  r.m.s.  value 


°26,70  '  1'689a"  ■ 


According  to  (5.18b),  we  should  have 


^2,25  +  a26,70^  ~  ^2,70 


(5.19a) 


This  relation  is  fulfilled  exactly,  since 


[(30. 2528m)2  +  (1.6898m)2:!15  =  30.3000m  . 


(5.19b) 


In  returning  to  the  theoretical  formulas  (5.17a,b),  except  for  we 
express  them  through  the  symbolism  similar  to  that  introduced  in  (5.18b): 


2 

an'+l,n" 


+  a 


n"+l ,  n 


=  a 


2 

n'  +1 ,  rf 


(5.20a) 


If  "n  -*•  <=°  ,  the  second  and  third  symbols  in  (5.20a)  become  what  is  called  the 
truncation  variance,  due  to  the  truncation  of  an  infinite  series  at  the  degrees 
n"  and  n',  respectively.  One  then  has 


a 


2 

n'+l,n" 


+  0 


n"  +  l,°° 


2 

a  ,  , 


(5.20b) 


5.3  Assessment  of  the  Numerical  Integration  Algorithm 


For  the  purpose  of  the  present  study,  two  equilateral  grids  are  generated, 
a  2°x2°  grid  and  a  4°x4°  grid.  In  referring  to  (5.4),  we  have  &4>  =  2°  for  the 
first  grid  and  A<j>  =  4°  for  the  second  grid;  the  intervals  AX  follow  from  the 

second  relation  in  (5.4).  Theoretically,  the  number  of  square  blocks  (P  ) 

2 

having  9  on  the  side  is  4tt/8  ;  upon  expressing  the  block  side  in  degrees,  this 
results  in 

Pt  =  (203. l°/0°) 2  .  (5.21) 

Accordingly,  the  theoretical  number  of  points  in  the  above  two  grids  is  2,578  and 
10,313,  respectively.  Due  to  the  rounding  of  the  number  of  grid  points  along 
each  parallel  to  the  nearest  integer  (for  safety  reasons  more  often  up  than  down), 
the  actual  total  number  (P)  of  the  grid  points  used  in  this  analysis  is  very 
slightly  higher.  The  complete  situation  is  summarized  by 


=  2°  ... 

, ..  P  =  10,313  , 

t 

P  *  10,328  ; 

(5.22a) 

-  4°  ... 

...  P  =  2,578  , 

P  =  2,586  . 

(5.22b) 

It  is  to  be  noted  that  at  both  poles  the  blocks  are  replaced  by  spherical  caps 
whose  size  is  about  25%  greater  here. 

Next,  we  turn  our  attention  to  the  number  of  coefficients  that  can  be  resolved 
from  a  9°x0°  equilateral  grid  of  geoid  undulations.  In  considering  that  e° 
represents  the  shortest  half-wavelength  of  the  geoidal  resolution  expressible  by 
an  (N,N)  set  of  S.H.  coefficients,  the  familiar  rule  of  thumb  suggests  that 

N  =  180°/ 9°  .  (5.23) 
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Since  the  total  number  of  coefficients  in  this  set  is  (N+l)  ,  in  conjunction 
with  the  above  two  grids  we  have 


8°  =  2°  .  (N.N)  =  (90,90)  ,  (N+l)z  =  8,281  ; 


(5.24a) 


0°  =  4°  .  (N.N)  =  (45,45)  ,  (N+l)^  =  2,116  . 


(5.24b) 


Upon  comparing  the  number  of  data  points  in  (5.22a,b)  with  the  number  of  unknown 
coefficients  in  (5.24a,b),  the  former  is  seen  to  exceed  the  latter  by  22-25%. 

Suppose  now  that  we  do  not  wish  to  accept  a  priori  the  limitations  (5.24a,b), 
and  attempt  to  establish  whether  a  larger  set  of  S.H.  coefficients  could  be 
resolved.  In  addressing  this  problem,  we  return  to  (5.7)  and  seek  to  solve  for 
the  S.H.  coefficients  as  unknown  parameters  in  a  least-squares  problem.  Since 
(5.7)  is  linear  in  these  coefficients,  an  observation  equation  can  be  formed 
with  zero  starting  values  for  AC  and  AS  : 

nm  nm 


v.  =  C  ...  RR  ( i )  ...  ;  ...  RS  (i)  ...  ] 

i  nm  nm 


-  Ni  * 


where  vi  is  the  residual,  Ni  is  the  (errorless)  geoid  undulation,  and  R  (i) 
and  S’  (i)  are  the  functions  from  (5.6c),  all  considered  at  the  point  i  with 

nm 

coordinates  (<j>i,Xi).  The  points  i,  i  =  1 ,  2,  ...  ,  P,  are  again  assumed  to  be 
distributed  in  an  equilateral  grid.  In  matrix  notations,  the  totality  of  these 


equations  is  represented  by 


V  =  A  X  -  CN .  3  . 

1 

Due  to  Ni  having  all  equal  weights,  the  regular  least-squares  solution  is 
X  =  (ATA)_1ATCN.]  . 

In  the  least-squares  (finite)  context,  the  consideration  of  crucial 
importance  is  again  whether  or  not  the  orthogonal  relations  (5.6a,b)  are 
reasonably  well  satisfied.  In  the  affirmative,  one  readily  confirms  that  the 

m 

matrix  of  normal  equations  (A  A)  is  in  theory  diagonal  with  all  the  diagonal 

elements  equal  to  R  P,  and  A  CN.3  is  a  vector  composed  of  two  distinct  groups 

of  elements,  the  first  group  having  the  form  R£N  R'  (i)  and  the  second  group 

i  nm 

having  the  form  RENjS  (i),  where  the  summation  £  applies  to  i  between  1  and 
P.  The  least-squares  solution  then  yields 


AC 

nm 

p 

R  (i) 

run 

AS 

nm 

=  ( l/R) ( 1/P)  £  N. 

i=l  1 

S  (i) 

nm 

which  is  precisely  (5.3).  And  in  the  negative,  ATA  can  no  longer  be  considered 
a  diagonal  matrix  and  the  least-squares  solution  no  longer  coincides  with  (5.3). 
It  would  now  involve  the  inversion  of  a  full  matrix  of  normal  equations,  and 
for  this  reason  alone  it  is  undesirable.  A  much  more  satisfactory  procedure  is 
to  simply  resort  to  a  finer  grid. 

As  we  have  seen  earlier,  a  coarse  grid  causes  the  coefficients  computed  >,y 
(5.3)  or  (5.5a,b)  to  be  distorted.  The  remedy  is,  of  course,  the  same  as  above- 
a  finer  grid.  We  can  conclude  this  discussion  by  stating  that  a  satisfactory 
solution  for  the  S.H.  coefficients  is  obtained  when  the  equilateral  grid  used 


P». 


in  their  evaluation  is  sufficiently  dense,  in  which  case  the  least-squares 
solution  coincides  in  practice  with  the  integrated  solution  (5.3)  or  (5.5a,b). 

A  denser  grid  leads  to  a  more  perfect  coincidence.  We  will  henceforth  consider 
only  the  integrated  solution  and  judge  its  quality  by  the  size  of  the  r.m.s. 
distortion  in  geoid  undulations  as  mentioned  earlier.  The  least-squares  analogy 
serves  merely  in  searching  for  the  largest  number  of  S.H.  coefficients  one  would 
ever  consider  resolving. 


Clearly,  if  the  number  of  observation  equations  were  smaller  than  the  number 
of  parameters,  a  unique  least-squares  solution  would  not  exist.  The  diagonal 
structure  of  the  matrix  of  normal  equations  would  be  destroyed  (the  latter  would 
be  singular)  and,  accordingly,  this  case  would  be  discarded.  This  means  that 
the  largest  number  of  parameters  we  shall  consider  is  equal  to  the  number  of 
grid  points.  If  the  latter  is  Pt,  and  the  largest  admissible  set  of  S.H.  coef¬ 
ficients  is  denoted  by  (N,  .  ,N.  .  ),  then  the  limiting  case  we  shall  still  examine 
is  characterized  by 


(N. .  +I)2  =  P  . 

lim  t 


In  conjunction  with  (5.21),  this  implies  that 


Nliin=  C203.1o/8°-i:int  .  (5.25) 

In  analogy  to  (5.24a,b)  one  can  now  write 


o 

CVJ 

II 

-  OWW  ■  <100'100>  • 

(Niim  +  D2  =  10,201  ; 

(5.26a) 
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o 

••  ■  < 49'  49>  • 

(Nlim*l)2  =  2,500. 

(5.26b) 

Upon  considering  P  from  (5.22a,b),  it  is  apparent  that  the  "redundancy"  has 
been  essentially  eliminated  (here  it  amounts  to  merely  1-3%  of  the  number  of 
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coefficients).  Ultimately,  computer  simulations  will  determine  whether  the 
number  of  coefficients  in  (5.26a,b),  or  at  least  that  in  (5.24a,b)  can  be 
satisfactorily  resolved  via  the  numerical  integration  (5.3)  or  (5.5a,b). 

The  computer  simulations  first  proceed  to  generate  errorless  Ni  according 
to  (5.1)  with  a  given  set  of  S.H.  coefficients  and  a  chosen  n',  then  to  compute 
new  S.H.  coefficients  within  (n'.n1)  via  (5.5a,b),  and  finally  to  recompute  fT 
with  the  latter  set  of  coefficients.  If  the  r.m.s.  of  the  distortions  in  N. 

i 

is  smaller  than  the  test  value  below,  the  algorithm  is  judged  satisfactory  with 
regard  to  the  (n',n‘)  set  of  coefficients  and  one  can  proceed  to  test  a  larger 
set.  The  same  procedure  applies,  of  course,  for  either  grid.  The  test  value 
(t)  is  determined  under  the  stipulation  that  the  r.m.s.  of  the  distortions  in 
Ni  should  increase  the  original  truncation  sigma  (c?n+1  J)  by  less  than  5%. 

This  occurs  for 

r.m.s.  (distortion)  <  t  -=  0.30  on,+1  ^  ;  (5.27a) 

in  particular, 


Co,  .  + 

n'  +1,°° 


(0.30)  o  ,  .  =  1.044  o  ,  .  , 

n'+l,00 


(5.27b) 


so  that  the  increase  in  the  truncation  sigma  is  actually  no  more  than  4.4%. 

In  order  to  find  t  needed  in  (5.27a),  the  truncation  sigma  is  computed  for 
various  n'.  A  practical  formula  giving  the  latter  can  be  transcribed  from  pages 
62  and  63  of  CBlaha,  19843  as  follows: 

oo 

V.i  - 5  D(0>  *  1  •  (5.28a) 

n  n=n' +1  n 

dn  =  0. 999617n+2  17,981  m2/C(n-l)(n-2) (n+24)]  . 


(5.28b) 


As  is  pointed  out  in  the  above  reference,  these  formulas  can  be  traced  to 
CTscherning  and  Rapp,  19743.  The  symbol  «  is  replaced  in  practice  by  a  suitable 
large  number,  such  as  1,000.  Equation  (5.28b)  presents  a  close  form  expression 
for  the  n-th  degree  variance  as  obtained  from  the  covariance  function.  A 
theoretical  formula  for  the  n-th  degree  variance  based  on  the  S.H.  potential 
coefficients  appeared  as  equation  (5.15).  The  values  of  t  computed  with  the  aid 


of  (5.28a,b)  are  listed  as  follows  (in  meters): 


a  m  4.89  3.03  2.00  1.80  1.66  1.63  1.20  1.06  0.94  0.85  0.77  0.71 

n'+l.oo' 

t  (m)  1.47  0.91  0.60  0.54  0.50  0.49  0.36  0.32  0.28  0.25  0.23  0.21 


The  actual  S.H.  potential  coefficients  used  in  the  simulations  of  geoid 
undulations  are  taken  from  a  (180,180)  set  computed  by  Rapp  at  The  Ohio  State 
University  (private  communication)  in  1980.  In  conjunction  with  the  2°x2° 
equilateral  grid,  the  r.m.s.  of  the  distortions  in  Ni  vary  according  to  selected 


degrees  n'  as  follows: 


r.m.s.  (m)  |  0.055  0.067  0.081  0.084  0.105  0.127  0.355  0.855 


With  regard  to  the  4°x4°  equilateral  grid,  the  simulations  yield 
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r.m.s.  m  0.113  0.126  0.129  0.180  0.247  0.555  2.389 


Distortions  very  similar  to  those  shown  in  (5.30)  and  (5.31)  have  been 
obtained  in  all  cases  where  the  testing  included  other  than  grid  points.  In 
one  example  of  this  kind,  geoid  undulations  were  tested  in  a  4.5°x4.5°  geo¬ 
graphical  grid,  where  the  integrated  solution  of  a  (49,49)  set  of  S.H.  coef¬ 
ficients  was  based  on  the  4°x4°  equilateral  grid.  Such  a  confirmation  of  the 
r.m.s.  of  the  distortions  is  reassuring,  especially  when  considering  that  the 
grid-point  distortions  in  (5.31)  are  already  deteriorating  for  this  truncation. 
Thus,  if  non-grid  points  should  exhibit  larger  distortions  than  the  grid  points, 
this  deterioration  would  have  certainly  manifested  itself  in  the  present 
example  chosen  for  its  increased  sensitivity. 

The  fact  that  the  distortions  beyond  n'  =100  in  (5.30)  and  beyond  n'  =49 
in  (5.31)  are  large  has  been  expected.  These  cases  entail  a  number  of  S.H. 
coefficients  larger  than  the  number  of  data  points,  i.e.,  exceed  the  limits 
imposed  by  (5.26a,b).  We  also  remark  that  when  one  attempts  to  resolve  a  larger 
set  of  coefficients  than  that  utilized  in  generating  the  errorless  data,  the 
r.m.s.  of  the  distortions  starts  a  slight  increasing  trend.  As  one  example,  an 
errorless  set  (50,50)  served  in  generating  Ni  in  the  2°x2°  equilateral  grid. 

The  r.m.s.  of  the  distortions  obtained  with  the  integrated  (50,50)  set  is  listed 
in  (5.30)  as  0.084  m.  But  when  the  numerical  integration  proceeded  through  a 
(60,60)  set,  this  r.m.s.  increased  to  0.094  m.  Although  the  integrated  values 
of  the  additional  coefficients  were  very  small  (ideally  they  should  be  zero), 
they  gave  rise  to  a  detectable  deterioration  in  the  results  nonetheless.  Such 
a  deterioration  appears  to  be  less  pronounced  with  the  2°x2°  than  with  the  4°*4° 
grid.  Far  from  being  worrisome,  it  is  consistent  with  the  trends  in  (5.30,31). 

In  both  the  2°x2°  and  4°x4°  equilateral  grids,  the  r.m.s.  of  the  distortions 
is  somewhat  smaller  when  computed  between  the  latitudes  +72°  than  when  computed 
for  the  whole  globe.  In  particular,  in  the  2°x2°  grid  it  is  reduced  to 


approximately  two-thirds  of  the  global  value  for  the  truncations  through  (90,90), 
and  in  the  4°*4°  grid  it  varies  approximately  from  one-half  to  two-thirds  of  the 
global  value  as  the  truncations  progress  from  (14,14)  to  (45,45).  This  improve¬ 
ment  is  attributed  to  a  more  regular,  nearly  square  shape  of  the  blocks  do. 
closer  to  the  equator,  as  opposed  to  the  blocks  near  the  poles  and  the  two 
spherical  caps  at  the  poles  themselves.  However,  this  approximately  one-third 
reduction  in  the  pertinent  r.m.s.  values  is  not  taken  into  account  in  the  analysis 
which  thus  might  be  slightly  conservative. 

Although  the  present  study  is  concerned  mainly  with  geoid  undulations, 
several  r.m.s.  values  of  the  distortions  in  gravity  anomalies  have  also  been 
computed.  In  analogy  to  (5.1)  and  (5.7),  including  the  same  assumptions,  error¬ 
less  values  of  Agi  are  first  generated  in  an  equilateral  grid  according  to  the 
formula 

4g  =  y  ?'  (n  -  1)  £  Rnm  +  4S,„  Sm)  , 

n=2  m =0 

where  y  is  the  average  value  of  gravity  (980  gal).  We  note  that  integral  formulas 
giving  the  S.H.  coefficients  based  on  gravity  anomalies  could  be  written  in 
analogy  to  (5.2)  and  (5.9),  except  that  y  and  Ag  would  replace  R  and  N,  respective 
ly,  and  l/(n-l)  would  be  factored  before  the  integral  sign.  The  same  changes 
would  take  place  also  with  regard  to  the  finite  case  (see  the  numerical  integra¬ 
tion  5.3  and  5.5a,b).  However,  this  modification  has  not  affected  the  present 
procedure  where  the  integrated  S.H.  coefficients  are  based  consistently  on  geoid 
undulations,  even  when  used  in  recomputing  Ag^  This  portion  of  the  simulations 
has  been  performed  only  in  the  2°*2°  equilateral  grid.  The  r.m.s.  distortions 
in  Agi  have  been  computed  in  conjunction  with  four  truncated  models  (n',n‘)  as 


follows: 


n' 


!  50  90  100  120 


(5.32) 


r.m.s.  (mgal)  i  0.3  1.1  4.8  13.6 

The  results  (5.29)  -  (5.32)  are  the  key  to  the  main  and  final  step  in  the 
analysis  —  the  determination  of  the  applicability  of  the  numerical  integration 
formulas  (5.3)  and  (5.5a,b).  This  determination  is  centered  on  geoid  undulations 
in  the  2°x2°  and  4°x4°  equilateral  grids,  and  thus  on  equations  (5.29)  -  (5.31) . 
Equation  (5.32)  serves  only  as  a  peripheral  confirmation  of  the  outcome  achieved 
with  the  2°x2°  grid.  A  graphic  representation  of  the  results  (5.29)  -  (5.31)  is 
offered  in  Fig.  1,  where  the  abscissa  displays  the  truncation  degrees  n',  and 
the  ordinate  depicts  either  the  values  of  t  from  (5.29)  or  the  r.m.s.  of  the 
distortions  in  geoid  undulations  from  (5.30)  and  (5.31).  The  former  case 
(equation  5.29)  is  identified  by  the  curve  bearing  the  description  "standard", 
and  the  latter  two  cases  (equations  5.30  and  5.31)  are  identified  by  the  curves 
bearing  the  descriptions  "2°x2°"  and  "4°x4°"  ,  respectively. 

If  we  were  to  proceed  strictly  according  to  the  philosophy  outlined  in  the 
paragraph  containing  (5.27a,b),  we  would  accept  the  truncation  (96,96)  for  the 
2°x2°  grid  and  the  truncation  (48,48)  for  the  4°x4°  grid  as  the  borderline  cases 
for  the  application  of  (5.3)  and  (5.5a,b).  These  truncations  would  also  approach 
the  limits  in  (5.26a,b),  respectively.  However,  Fig.  1  reveals  an  unusual  aspect 
conveying  an  important  message  in  this  respect.  Both  curves  "2°x20"  and  "4°x4°" 
indicate  in  a  clearcut  fashion  that  in  this  analysis,  the  rule  embodied  by 
(5.24a,b)  is  much  more  significant  than  might  be  expected  from  an  approximate 
rule  of  thumb. 

Indeed,  in  the  important  case  "2°x2°"  the  r.m.s.  of  the  distortions  in  geoid 
undulations  through  the  truncation  (90,90)  remains  comfortably  low,  well  below 
the  standard  t  irs  Fig.  1.  Within  this  range,  the  r.m.s.  increases  essentially 
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Fig.  1 

Display  of  the  curve  representing  the  standard  of  comparison  (t)  for 
different  truncation  degrees  n',  and  the  curves  representing  the  r.m.s 
of  the  distortions  in  geoid  undulations  due  to  numerical  integration 
in  2°x2°  and  4°x4°  equilateral  grids  for  different  n' 


in  a  linear  fashion  and  at  a  very  low  rate.  However,  beyond  the  truncation 
(90,90)  the  r.m.s.  deteriorates  rapidly,  as  the  figure  displays  in  a  self-evident 
manner.  The  limited  analysis  of  the  distortions  in  Ag  summarized  in  (5.32) 
corroborates  this  finding.  In  particular,  if  the  r.m.s.  values  from  (5.32)  were 
included  in  the  figure,  they  would  produce  a  pattern  almost  identical  to  that  of 
the  curve  "2°x2°"  (except,  perhaps,  that  the  slope  of  the  deterioration  beyond 
n'=90  would  be  slightly  lower).  The  curve  depicted  as  "4°x4°"  in  Fig.  1  also 
has  clearcut  characteristics,  very  similar,  in  fact,  to  those  of  the  "2°x20" 
curve.  The  truncation  beyond  which  a  sharp  deterioration  occurs  is  now  (45,45), 
again  demonstratively  in  line  with  the  rule  of  thumb. 

In  adopting  n'=90  in  the  2°x2°  case  and  n'=45  in  the  4°x4°  case,  the  r.m.s. 
of  the  distortions  become  0.127  m  and  0.247  m,  respectively  (see  equations  5.30 
and  5.31).  These  values  amount  to  merely  0.14an,+1  ^  in  both  cases  (see 
equation  5.29).  Accordingly,  the  increase  in  the  original  truncation  sigma  is 
not  4.4%  as  in  (5.27b),  but  only  1%,  since 


Ca2  lx.  +  (0.14  a  ,  ,  )V  »  1.01  o 

n'+l,®  1  n ' +1,»/  n  +1, 1 


(5.33) 


We  may  thus  conclude  that  only  negligibly  small  worsening  (1%)  in  the  original 
truncation  sigma  occurs  due  to  the  numerical  integration  (5.3)  and  (5.5a,b), 
provided  the  truncated  model  is  (90,90)  for  the  2°x2°  equilateral  grid,  and 
(45,45)  for  the  4°x4°  equilateral  grid. 

This  conclusion  embodies  two  encouraging  signs.  The  first  points  to  the 
truncation  degrees  as  being  comfortably  high,  and  the  second,  already  mentioned, 
points  to  these  degrees  as  corresponding  exactly  to  what  has  been  termed  the  rule 
of  thumb  in  relating  the  gravity  field  resolution  to  the  size  of  a  set  of  S.H. 
potential  coefficients  used  in  the  field's  representation.  Although  this  out¬ 
come  has  been  produced  with  computer  simulations  performed  in  the  spherical 
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approximation,  detailed  determinations  of  the  actual  gravity  field  can  follow 
the  same  principles  and  take  advantage  of  the  same  computational  algorithms. 
Since,  at  AF6L,  the  geoid  determination  from  satellite  altimetry  is  based  on  a 
(14,14)  reference  surface  rather  than  an  ellipsoid,  the  spherical  approximation 
entails  errors  on  the  order  of  a  few  centimeters  (instead  of  decimeters),  and  is 


thus  deemed  inconsequential. 


5.4  Smoothing  Effect 


We  have  seen  that  the  numerical  integration  (5.5a,b)  leads  to  a  successful 
modeling  of  geoid  undulations,  and  thus  of  the  disturbing  potential,  provided 
it  is  used  up  to  the  degree  and  order  (N,N),  where  N=18OC)00  according  to  the 
familiar  rule  of  thumb  presented  as  (5.23),  in  which  6°  represents  the  side  of 
equilateral  blocks  containing  data  points.  Except  for  regions  near  the  poles, 
the  blocks  are  essentially  square  and  the  data  points  are  located  at  their 
centers.  The  latter  are  conceived  as  forming  an  equilateral  grid  according  to 
(5.4),  where  A<J>  corresponds  to  0  above. 

The  distortions  in  geoid  undulations  (differences  between  the  errorless  and 
the  modeled  values)  have  been  judged  in  the  context  of  one  and  the  same  truncat¬ 
ed  model,  i.e.,  in  the  situation  where  both  the  errorless  and  the  modeled  values 
are  computed  within  the  same  degree  and  order  truncation  (n',n‘).  However,  if 
the  errorless  values  are  generated  with  an  "original"  set  (N,N)  of  spherical- 
harmonic  (S.H.)  coefficients,  and  the  modeled  values  are  computed  with  an  inte¬ 
grated  set  (n'.n1),  where  n'<N,  the  distortions  owe  their  existence  primarily 
to  a  smoothing  effect.  One  can  then  judge  the  quality  of  the  modeled,  smoothed 
geoid  (n',n')  by  comparing  it  with  an  errorless  geoid  (n',n')  generated  with  the 
(n',n')  subset  of  the  original  set  of  S.H.  coefficients. 

Ideally,  the  above  two  surfaces  should  coincide,  i.e.,  the  geoidal  detail 
due  to  the  degrees  n'+l  through  N  should  be  suppressed.  This  desirable 
characteristic  hinges  on  the  fulfillment  of  the  orthogonal  reations,  as  can  be 
shown  in  theory  with  the  aid  of  (5.13).  If  we  subject  this  equation  to  the 
procedure  which  led  from  (5.7)  to  (5.8),  its  right-hand  side  becomes  zero  for 
any  s<_n'  due  to  (5.6a).  But  this  implies  that 
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//N(n")  _sr  do  =  JjN( n ' )  _sr  da  , 


r  <  s  <  n'  <  n"  . 


(5.34) 


Therefore,  when  the  S.H.  coefficients  are  evaluated  through  the  degree  and  order 
(n'.n1)  as  in  (5.9),  N  under  the  integral  sign  can  be  expressible  within  any 
model  (n",n")  as  long  as  n">n'.  In  other  words,  the  geoidal  detail  described 
by  N(n")-N(n‘)  has  no  bearing  on  the  theoretical  values  of  the  S.H.  coefficients 
within  the  truncation  (n',n'). 

Computer  simulations  have  confirmed  the  above  outcome  in  the  finite  case. 

The  results  presented  below  are  sufficient  to  illustrate  this  fact.  They  have 
been  obtained  with  the  aid  of  a  2°x2°  equilateral  grid  of  geoid  undulations 
generated  via  S.H.  coefficients  complete  through  the  degree  and  order 
(N,N)=(90,90) ,  where  N  satisfies  the  rule  (5.23).  Of  the  models  (n',n'),  where 
the  smoothing  effect  takes  place  for  any  n'<N,  three  are  selected  for  a  limited 
analysis,  namely  (25,25),  (45,45),  and  (70,70). 

In  each  of  the  three  cases  (n',n'),  three  values  r.m.s. ^  i=l,2,3,  are  listed: 


n’ 

_ 1 

25 

45 

70 

r.m.s. 

(m) 

0.067 

0.081 

0.105 

r.m.s. 2 

(m) 

0.068 

0.082 

0.105 

r.m.s. 3 

(m) 

0.002 

0.002 

0.003 

(5.35) 


The  first  row  of  results,  i=l,  contains  the  entries  from  the  corresponding  places 
in  (5.30)  and  is  listed  merely  for  the  sake  of  comparison.  The  entries  depict 
the  r.m.s.  of  the  discrepancies  in  geoid  undulations  between  the  (n',n')  integrat¬ 
ed  values  (i.e.,  geoid  undulations  computed  with  the  n' ,n'  set  of  S.H.  coefficients 
evaluated  via  the  numerical  integration  5.5a,b)  based  on  the  (n',n')  errorless 
values  (i.e.,  geoid  undulations  N  in  5.5a)  on  one  hand,  and  the  (n',n‘)  errorless 
values  themselves  on  the  other  hand. 


The  second  row  of  results,  i=2,  contains  the  r.m.s.  of  the  discrepancies 
between  the  (n',n')  integrated  values  based  on  the  (90,90)  errorless  values, 
and  the  (n',n‘)  errorless  values.  The  (n',n')  integrated  values  in  this  case 
encompass  the  smoothing  effect.  As  is  apparent  from  this  row,  the  smoothed 
{ n '  , n ' )  values  approximate  the  errorless  values  as  closely  as  the  integrated 
values  based  directly  on  the  smooth  geoid  (see  the  case  i=l).  This  is  strongly 
evidenced  by  the  last  row,  i=3,  depicting  the  r.m.s.  of  the  differences  between 
two  kinds  of  (n’.n1)  integrated  values,  the  first  based  on  the  (n',n‘)  errorless 
geoid  and  the  second  based  on  the  (90,90)  errorless  geoid. 

In  complete  analogy  to  (5.35),  the  following  results  are  obtained  with  the 
aid  of  a  4-4"  equilateral  grid  of  geoid  undulations  generated  via  an 
(N  ,N)  =  (45 ,45)  set  of  S.H.  coefficients  where  the  models  (n',n')  are  selected  as 
(14,14),  (25,25),  and  (36,36): 


n' 

14 

25 

36 

r.m. s. 

(m) 

0.113 

0.129 

0.180 

r.m.s. 2 

(m)  j 

0.113 

0.130 

0.181 

r.m.s. 

(m) 

0.004 

0.005 

0.010 

The  row  i=l  contains  the  corresponding  entries  from  (5.31),  pertaining  to  the 
differences  between  the  ( n 1  ,n ' )  integrated  values  based  on  the  (n',n')  errorless 
values,  and  the  ( n '  ,n ' )  errorless  values  themselves.  The  row  i=2  pertains  to  the 
differences  between  the  (n',n')  integrated  values  based  on  the  (45,45)  errorless 
values,  and  the  (n',n')  errorless  values.  And  the  row  i=3  pertains  to  the  dif¬ 
ferences  between  the  above  two  kinds  of  (n',n‘)  integrated  values. 


The  results  in  (5.36)  follow  a  very  similar  pattern  to  those  in  (5.35). 

In  comparing  the  corresponding  cases  (25,25),  one  notices  that  all  three 
entries  in  (5.36)  are  approximately  the  double  of  their  counterparts  in  (5.35). 
This  tendency  can  be  observed  already  in  (5.31)  versus  (5.30)  for  relatively 
low  truncation  degrees  n' ,  including  the  truncation  (30,30)  for  which  the 
respective  r.m.s.  values  would  be  approximately  0.071  m  and  0.152  m.  Clearly, 
this  one-half  reduction  in  the  discrepancy  values  is  imputable  to  the  two-fold 
increase  in  density  of  the  equilateral  grid  along  each  dimension. 

One  important  conclusion  of  this  analysis  consists  in  the  finding  that  the 
portion  (5.6a)  of  the  orthogonal  relations  is  satisfied  particularly  well  with 
the  selected  truncation  degrees  n'<N,  quite  representative  of  realistic  levels 
of  smoothing.  Indeed,  the  small  magnitude  of  the  values  in  the  last  rows  of 
(5.35)  and  (5.36)  indicates  that  the  effect  of  geoid  undulations  due  to  the 
degrees  n'+l  through  N  is  almost  perfectly  suppressed  in  the  application  of  the 
numerical  integration  algorithm. 


6.  APPLICATIONS  OF  THE  NUMERICAL  INTEGRATION  ALGORITHM 
6 . 1  Global  Representation  of  the  Gravity  Field 

A  high  degree  and  order  S.H.  expansion  of  the  disturbing  potential  provides 
for  a  detailed  representation  of  the  earth's  gravity  field  and  its  fundamental 
surface,  the  geoid.  If  the  reference  field  is  defined  through  a  S.H.  expansion 
of  degree  and  order  substantially  higher  than  two  (which  corresponds  essentially 
to  an  ellipsoidal  field),  the  spherical  approximation  is  acceptable.  Thus,  when 
referring  to  this  field,  geoid  undulations  and  other  functions  of  the  geopoten¬ 
tial  (gravity  anomalies,  deflections  of  the  vertical,  etc.)  can  be  regarded  as 
functions  on  a  sphere  and  expressed  in  terms  of  surface  spherical  harmonics. 

This  has  been  the  case  at  AFGL,  where  the  global  first-phase  altimeter  adjustment 
has  determined  a  (14,14)  reference  field  as  a  basis  for  subsequent  representa¬ 
tions  of  the  gravity  field  detail,  especially  the  oceanic  geoid. 

As  we  have  seen  in  Chapter  5,  if  geoid  undulations  are  supplied  in  a  suf¬ 
ficiently  dense  equilateral  grid  all  o'^er  the  globe,  the  numerical  integration 
algorithm  permits  a  direct  evaluation  of  the  desired  S.H.  potential  coefficients 
in  the  context  of  the  above  spherical  approximation.  This  set  can  serve  in 
predicting  the  desired  geophysical  quantities  (especially  the  geoid  undulations 
themselves)  at  any  location.  The  prediction  grid  needed  for  the  construction 
of  contour  maps  can  then  be  made  geographic  and  can  be  densified  at  will. 

One  major  source  of  geoid  undulations  distributed  in  an  equilateral  grid 
over  the  world's  oceans  is  the  modified  collocation  with  noise,  whose  main 
features  were  recapitulated  in  Chapter  4.  The  grid  density  determines  the 
largest  size  (N,N)  of  the  set  of  S.H.  coefficients  which  can  be  satisfactorily 
resolved  via  numerical  integration  (5.5a,b).  This  relationship  was  demonstrated 


in  Section  5.3,  resulting  in  the  confirmation  of  the  familiar  rule  of  thumb 
(5.23).  Thus,  a  2°x2°  equilateral  grid  of  geoid  undulations  allows  for  a 
reliable  evaluation  of  a  (90,90)  set  of  S.H.  coefficients.  The  same  grid  can 
be  used  in  expressing  a  smoother  version  of  the  gravity  field  simply  by  trunca¬ 
ting  the  (90,90)  set  at  a  lower  degree  and  order. 

In  applying  the  numerical  integration  algorithm,  geoidal  values  over  the 
continents  must  also  be  given.  The  most  expedient  procedure  is  to  simulate  them 
using  some  a  priori  set  of  S.H.  potential  coefficients  complete  through  the 
required  degree  and  order  (e.g.,  a  90,90  set  as  above).  This  set  need  not  be 
of  high  accuracy,  since  the  most  important  product  of  satellite  altimetry  is  the 
determination  of  the  oceanic  geoid;  the  contours  over  the  land  masses  can  ulti¬ 
mately  be  disregarded.  On  the  other  hand,  the  simulated  values  must  be  realistic 
(including  the  appropriate  geoidal  roughness)  because  otherwise  the  oceanic  con¬ 
tours  within  some  distance  from  the  coasts  would  be  deformed.  We  can  conclude 
by  stating  that  the  oceanic  geoid  expressed  by  the  new,  integrated  set  of  S.H. 
coefficients  represents  an  improvement  over  that  expressed  by  the  a  priori  set 
not  benefiting  from  the  contribution  of  altimeter  data. 


6.2  Global  Representation  of  Tidal  Effects 


A  global  tidal  adjustment  within  the  first-phase  adjustment  of  satellite 
altimetry  (including  S.H.  potential  coefficients  as  parameters)  was  described 
in  CBlaha,  1982]  and  Chapter  3  of  CBlaha,  1984D.  Subsequently,  an  optional 
tidal  (re)adjustment  within  the  second-phase,  large-area  adjustment  of  satellite 

altimetry  (including  point-mass  magnitudes  as  parameters)  was  described  in 
Chapter  5  of  CBlaha,  1984].  Since,  in  the  second  phase,  the  coefficients  of 
tidal  parameters  in  the  observation  equations  are  adopted  from  the  first  phase, 
the  second-phase  adjustment  will  no  longer  be  mentioned  in  this  analysis. 

The  a  priori  tidal  information  in  the  first-phase  adjustment  at  AFGL  has 
been  supplied  via  S.H.  tidal  coefficients.  Since  the  Legendre  polynomials  and 
the  associated  Legendre  functions  as  well  as  the  multiple-angle  trigonometric 
functions  in  this  adjustment  are  computed  within  an  (n',n')  S.H.  model  even 
without  tidal  considerations ,  little  additional  effort  is  needed  when  the  a  priori 
tidal  effects  are  expressed  by  a  S.H.  expansion  within  the  same,  or  lower,  degree 
and  order  truncation.  In  recent  adjustments  of  satellite  altimetry,  the  S.H. 
expansion  in  the  first  phase  has  been  performed  in  a  (14,14)  truncated  model. 

On  the  other  hand,  the  information  pertaining  to  seven  diurnal  and  semidiurnal 
tidal  constituents  has  been  included  via  S.H.  tidal  coefficients  in  two  (12,12) 
sets  per  constituent,  as  is  described  in  CBlaha,  1982]  in  conjunction  with 
CEstes,  1980]. 

However,  much  more  complete  (in  terms  of  a  dense  global  coverage)  and  more 
accurate  tidal  information  has  become  recently  available  in  various  NSWC  Techni¬ 
cal  Reports  by  E.  W.  Schwiderski.  The  relative  accuracy  for  each  ocean-tide 
constituent  treated  is  stated  to  be  5  cm  in  the  open  oceans.  The  tidal  amplitudes 
and  Greenwich  phase  angles  of  these  constituents  are  tabulated  in  a  1°*1° 
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geographical  grid  for  all  the  oceanic  areas  in  the  world,  and  this  information 
is  also  available  in  the  form  of  global  corange  and  cotidal  maps.  Such  know¬ 
ledge  allows  one  to  compute  the  ocean  tide  as  a  function  of  time  and  location, 
constituting  a  stepping  stone  in  the  computation  of  the  surface  tide  directly 
related  to  altimeter  measurements  (see,  e.g..  Chapter  3  of  CBlaha,  19843). 


The  ocean-tide  information  for  each  tidal  constituent  can  be  obtained  from 
the  NSWC  reports  as  just  mentioned,  or,  less  directly,  from  the  S.H.  tidal  coef¬ 
ficients  generated  in  turn  from  the  information  contained  in  these  reports. 

Since  the  present  procedure  at  AFGL  allows  for  an  actual  adjustment  of  tidal 
amplitudes  and  Greenwich  phase  angles  as  parameters,  provided  the  tidal  informa¬ 
tion  is  supplied  via  S.H.  tidal  coefficients,  the  latter  approach  offers  certain 
advantages.  In  particular,  it  permits  an  analysis  of  how  well  the  tidal  informa¬ 
tion  is  consistent  with  satellite  altimeter  observations  (one  could  judge  by  the 
corrections  attributed  to  the  tidal  parameters  by  the  adjustment  versus  the 
a  priori  sigmas,  etc.).  Such  an  analysis  can  indicate  possible  weaknesses  in 
the  tidal  information,  the  presence  of  systematic  orbital  errors  of  tidal 
frequencies,  etc. 

The  algorithm  for  computing  the  ocean  tide  based  on  the  information  contained 
in  the  NSWC  reports  is  presented  as  follows.  First,  we  recapitulate  the  basic 
formulas  pertaining  to  a  constituent  "j"  from  pages  48  and  49  in  CBlaha,  19823, 
omitting  the  subscript  j: 


«  =  +  ^  • 


(6.1a) 


C,  =  cosa  A  cos^  ,  £ „  =  sina  A  sin^  , 


(6.1b) 


£  =  £(4>,A,t)  =  ocean  tide, 
a  =  a(t)  =  Greenwich  argument, 

A  =  A( 4> , x)  =  tidal  amplitude, 

1 1>  =  ij>(4>»A)  =  Greenwich  phase  angle, 

and  where 

4>,A  =  geocentric  latitude  and  longitude, 
t  =  time. 

Next,  we  form  F  and  G,  regarded  as  functions  on  a  sphere,  such  that 

F  =  F(<M)  =  A(4>,A)  cosij;(<M)  . 

G  =  G(4>,A)  =  A(4>,A)  sin4/(<f>,A)  . 

In  expressing  these  functions  in  terms  of  surface  spherical  harmonics,  we  have 

°°  n  _  _  _ 

F  h  A  cos<p  =  Z  E  (a  cos  mA  +  b  sin  mA)P  (sin4>)  ,  (6.2 

_  „  run  nm  nm 

n=0  m=0 

00  n  _  _  __ 

G  e  A  s i nip  =  Z  Z  (c  cos  mA  +  d  sin  mA)P  (si n<j> )  .  (6.2 

.  _  nm  nm  nm 

n=0  m=0 

The  tidal  amplitudes  A  (in  cm)  and  Greenwich  phase  angles  \p  (in  deg.)  are 
available  from  the  NSWC  reports.  For  example,  CSchwiderski ,  1979D  presents 
these  values  for  the  semidiurnal  principal  lunar  tide  M2-  One  can  thus  form  F 
and  G  on  the  left-hand  sides  of  (6.2a,b)  for  any  oceanic  point  by  finding  its 
A  on  the  corange  map  (or  in  the  appropriate  table  of  ocean  tide  amplitudes) 
and  its  ip  on  the  cotidal  map  (or  in  the  corresponding  table  of  ocean  tide 
Greenwich  phases). 


In  order  to  obtain  the  S.H.  tidal  coefficients  a  ,  E  ,  c  ,  and  3 

nm  run  nm  run 

via  integral  formulas  presented  below,  the  values  of  F  and  6  must  be  known 
(in  an  equilateral  grid)  over  the  whole  globe,  not  just  over  the  oceanic  areas. 
Similar  to  the  previous  section,  values  associated  with  the  land  masses  must  be 
supplied  in  some  manner.  Stipulating  unrealistic  values,  such  as  zeros,  would 
be  harmful  to  the  tidal  representation  in  the  coastal  regions. 

This  problem  can  be  solved  as  follows.  First,  similar  to  the  approach  of 
Chapter  5  (treating  N  as  point  values)  we  compute  the  point  values  of  the 
functions  F  and  G  over  the  oceanic  areas  in  a  suitable  equilateral  grid  consistent 
with  the  desired  degree  and  order  truncation  (n',n‘).  As  an  example,  a  4°x-4° 
equilateral  grid  might  be  sufficiently  dense  because  these  coefficients  would 
seldom  be  needed  beyond  the  (45,45)  truncation.  Lowering  the  degree  and  order 
of  the  truncation  from  this  level  would  amount  to  further  smoothing  of  the  infor¬ 
mation  contained  in  the  NSWC  reports. 

Next,  we  draw  contours  depicting  the  functions  F  and  G,  completing  them 
visually  over  the  land  masses  in  a  manner  compatible  with  the  oceanic  contours. 

The  land  values  of  F  and  G  have  no  physical  meaning,  and  no  use  of  their  own, 
other  than  to  make  possible  an  evaluation  of  S.H.  tidal  coefficients.  The  latter 
will  serve  to  reproduce  the  oceanic  values  and  will  never  be  used  in  conjunction 
with  any  point  over  land. 

Finally,  having  the  functions  F  and  G  in  (6.2a,b)  available  in  an  equilateral 
grid,  we  proceed  to  the  evaluation  of  the  desired  S.H.  tidal  coefficients  in 
analogy  to  (5.5a,b) : 
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(6.3a) 


(6.3b) 


(6.3c) 


where  K  is  again  the  number  of  parallels  containing  grid  points  and  pk  is  the 
number  of  grid  points  on  the  parallel  "k".  Since  (6.3b)  differs  from  (6.3a) 
merely  in  the  value  G  replacing  F  for  each  point  of  the  (common)  grid,  these  two 
relations  are  evaluated  simultaneously. 


The  resulting  set  (n',n')  of  S.H.  tidal  coefficients  can  now  be  used  in 
(6.2a,b),  where  the  symbol  «°  is  replaced  by  n' ,  in  order  to  recompute  the 
functions  F  and  G  at  the  grid  points  (for  verification  purposes),  as  well  as  to 
evaluate  them  at  any  other  locations.  Most  important  in  this  respect  are  points 
in  the  oceanic  regions  containing  altimeter  data.  In  utilizing  (6.2a,b)  in 
(6.1a,b),  the  ocean  tide  for  the  pertinent  constituent  can  be  expressed  as 


=  COSu 


+  sina 


n '  n 
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(a  cos  mx  +  b  sin  mX)P  ( sin<t>) 
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(6.4) 


Such  values  of  the  ocean  tide  for  the  constituent  "j",  improved  by  comparison  to 


the  previous  procedure,  can  serve  in  the  tidal  adjustment  exactly  as  £  did  in 
CBlaha,  19823  in  conjunction  with  Chapter  3  of  CBlaha,  19843. 


7.  CONCLUSIONS 


Two  second-phase  techniques  have  been  developed  in  recent  years  at  AFGL 
based  on  the  results  of  the  global  spherical -harmonic  (S.H.)  treatment  of  alti- 
metric  data.  One  of  these  techniques  has  been  documented  as  the  point-mass 
(P.M.)  adjustment,  and  the  other  has  been  described  under  the  name  of  modified 
collocation  with  noise.  The  former  is  characterized  by  a  simultaneous  adjustment 
of  all  the  P.M.  parameters;  however,  the  matrix  of  normal  equations  to  be  in¬ 
verted  can  be  made  banded  or,  in  the  presence  of  (optional)  tidal  parameters, 
banded-bordered.  And  the  latter  is  noted  for  its  economical  property  of  avoiding 
an  inversion  of  a  large  matrix  (such  as  is  performed  in  the  P.M.  adjustment), 
proceeding  instead  to  the  inversion  of  a  number  of  relatively  small  matrices,  one 
per  observation  point;  there  is  no  limit  to  the  number  of  these  points.  The  term 
"modified"  attributed  to  the  collocation  technique  does  not  concern  the  philos¬ 
ophy  of  the  least-squares  collocation  with  noise  per  se,  but,  rather,  its  specific 
application  aimed  at  describing  a  smoothed-out  gravity  field,  in  which  the  part  of 
the  signal  beyond  the  desired  smoothing  level  is  pushed  into  the  realm  of  "noise". 

The  P.M.  adjustment  has  been  designed  in  two  versions,  the  single-layer  mode, 
where  all  the  point  masses  are  situated  at  one  depth  below  the  surface  of  the 
reference  ellipsoid,  and  the  double-layer  mode,  where  to  each  such  point  mass  is 
attributed  another  one  below  it,  separated  from  it  by  the  same  vertical  distance 
v.  The  number  of  parameters  in  the  latter  adjustment  setup  is  unchanged  from  the 
former,  due  to  the  stipulation  that  the  twin  point  masses  differ  only  in  sign, 
not  in  magnitude.  Even  so,  the  double-layer  mode  is  more  time-consuming  because 
of  an  increased  number  of  computer  operations  needed  in  the  formation  of  observa¬ 
tion  equations. 


From  the  tests  conducted  in  Chapter  3  it  transpires  that  in  some  cases 
the  double-layer  adjustment  may  accommodate  the  observations  (taken  as  minus 
the  geoidal  residuals  from  the  global  S.H.  adjustment)  slightly  better  than  its 
single-layer  counterpart.  However,  the  differences  are  small  and  the  improve¬ 
ments  (in  terms  of  the  r.m.s.  values)  do  not  exceed  7%  for  any  choice  of  v. 

Since  variations  in  the  depth  of  the  shallower  layer  itself  are  also  shown  to 
have  only  a  small  effect  on  the  resulting  r.m.s.  of  the  residuals,  it  follows 
that  the  power  of  the  resolution  is  linked  almost  entirely  to  the  horizontal  dis¬ 
tribution  of  the  point  masses.  It  can  be  concluded  that  due  to  its  economic 
superiority,  the  single-layer  mode  is  preferable  to  its  double-layer  counterpart, 
especially  if  the  P.M.  adjustment  should  be  performed  on  the  global  oceanic  scale. 

If  a  very  detailed  gravity  field  resolution  is  needed  in  areas  of  special 
interest,  one  is  compelled  to  consider  an  additional  treatment  of  altimeter  data, 
such  as  a  third-phase  P.M.  adjustment.  An  adjustment  of  this  type  is  the  topic 
of  Chapter  2.  It  is  based  on  the  second-phase  results  in  much  the  same  way  as  the 
second-phase  P.M.  approach  is  based  on  the  outcome  of  the  first-phase  (S.H.) 
adjustment.  In  the  third  phase,  the  point  masses  in  an  individual  area  of  interest 
form  a  substantially  denser,  and  commensurably  shallower  grid  than  that  character¬ 
izing  the  underlying  second-phase  P.M.  adjustment.  The  final  quantities  are  ob¬ 
tained  as  the  sum  of  three  parts,  each  corresponding  to  the  appropriate  adjustment 
phase.  These  quantities  can  subsequently  be  used  in  the  construction  of  contour 
maps  representing  the  geoid  and  other  functions  of  the  geopotential  (such  as 
gravity  anomalies  or  deflections  of  the  vertical).  In  basing  one  adjustment 
phase  on  the  residuals  from  the  previous  phase,  the  concept  of  point  masses 
ensures  a  gradual  transition  of  contour  lines  between  the  lower  and  the  higher 
resolution  areas. 


The  same  chapter  also  addresses  a  third-phase  solution  based  on  the 
second-phase  results  in  terms  of  the  modified  collocation  with  noise.  Except 
for  changes  at  the  algorithmic  level,  such  a  procedure  is  governed  by  the  same 
principles  as  the  third-phase  versus  the  second-phase  P.M.  approach.  As  an 
example  of  the  algorithmic  changes,  the  (n,n)  S.H.  model  used  in  the  second 
phase  as  "reference  field"  is  replaced  by  an  (n*  »n ' )  model  in  the  third  phase, 
which  is  precisely  the  model  obtained  as  a  solution  of  the  second  phase. 
Accordingly,  the  covariance  and  cross-covariance  functions  in  the  third  phase 
exhibit  sharper,  more  localized  characteristics  compared  to  their  counterparts 
in  the  second  phase. 

Since  the  collocation  approach  does  not  provide  for  adjustment  of  tidal 
parameters,  the  third-phase  collocation  solution  could  be  combined  with  the 
second-phase  P.M.  adjustment  including  tidal  effects.  In  this  way,  advantages 
of  both  methods  would  be  allowed  to  manifest  themselves.  The  second  phase, 
carried  out  with  a  relatively  sparse  grid  of  point  masses,  would  offer  the  benefit 
of  tidal  adjustment  at  the  level  of  individual  ocean  basins,  while  the  third 
phase,  performed  in  the  collocation  mode,  would  produce  a  detailed  resolution  of 
the  gravity  field  on  a  desired  scale,  an  ocean-basin  or  even  a  global  scale. 

Due  to  a  similarity  between  the  surface  spherical-harmonic  expansion  of  the 
ocean  tide  for  a  constituent  "j"  and  the  expansion  of  the  geoid  undulation,  one 
could  express  the  former  in  terms  of  "tidal"  point  masses,  just  as  the  latter 
has  been  expressed  in  terms  of  "regular"  point  masses.  This  possibility  is  ex¬ 
plored  in  Appendix  2.  A  practical  procedure  which  would  treat  the  new  P.M. 
parameters  together  with  the  existing  ones  faces  one  substantial  difficulty  -- 
the  computer  core  limitations.  In  particular,  if  the  present  P.M.  adjustment 
should  include  tidal  point  masses,  their  number  would  have  to  be  small  and, 
therefore,  their  distribution  sparse,  to  the  point  of  being  equivalent  to  a 


(22,22)  S.H.  expansion.  Such  a  resolution  power  is  only  slightly  better  than 
that  provided  by  the  existing  capability  of  the  first-phase  adjustment  alone, 
although  the  latter  affects  only  two  parameters  per  tidal  constituent  (the  tidal 
amplitude  and  the  Greenwich  phase  angle).  However,  these  parameters  have  a 
clearcut  physical  meaning  and  a  number  of  advantages,  as  is  discussed  in  Section 
2.1  of  CBlaha,  1983].  Be  that  as  it  may,  the  resolution  power  per  se  is  not 
affected  by  a  particular  mode  of  tidal  adjustment. 

The  present  second-phase  P.M.  adjustment  allows  for  a  supplementary  tidal 
adjustment  of  the  same  two  parameters  per  constituent  as  discussed  above.  The 
adjustment  corrections  apply  to  the  tidal  effects  in  predetermined  regions, 
typically  in  individual  ocean  basins,  and,  similar  to  the  first  phase,  lend 
themselves  readily  to  a  physical  interpretation.  It  is  especially  noteworthy 
that  the  three  tidal  constituents  M2,  N2,  and  0A  considered  in  view  of  SEASAT 
altimetry  add  only  a  negligible  computational  burden  to  the  adjustment  process 
by  bringing  six  parameters  into  the  border  of  normal  equations;  their  solution 
is  carried  out  efficiently  by  the  algorithm  presented  in  Chapter  5  of  CBlaha 
1984],  This  property  is  important  when  the  computer  core  is  already  utilized 
to  its  capacity,  and  leads  to  the  conclusion  that  the  computer  limitations  make 
it  impractical,  at  the  present,  to  replace  the  supplementary  tidal  adjustment 
within  the  second-phase  P.M.  approach  by  a  new  procedure  introducing  a  number 
of  tidal  P.M.  parameters  into  the  adjustment.  However,  an  algorithm  to  this 
effect  is  developed  in  Appendix  2,  which  can  be  put  to  use  when  permitted  by 
computational  and  economic  factors. 

In  the  past  treatment  of  satellite  altimetry,  the  first-  and  even  the  second- 
phase  adjustment  exhibited  .large  negative  altimeter  residuals  in  areas  where 
satellite  arcs  cross  an  important  trench,  such  as  the  Puerto  Rico  trench.  Thus, 
altimeter  residuals  may  have  a  significant  role  to  play  in  the  detection  of 
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bathymetric  and  other  oceanic  anomalies  (e.g.,  geomagnetic  anomalies).  This 
topic  is  discussed  at  length  in  Appendix  3.  Since  it  appears  that  a  second- 
phase  approach  would  be  insufficient  for  a  reliable  detection  of  even  some  im¬ 
portant  trenches,  and  a  third-phase  approach  would  represent  a  lengthy  and 
uneconomical  undertaking  for  this  purpose,  one  is  well  advised  to  turn  to  the 
first-phase  residuals  as  a  basic  tool  of  such  an  endeavor.  However,  instead 
of  considering  only  the  actual  residuals  at  observation  points  selected  at  - 
or  ^-intervals  in  the  first-phase  adjustment  of  SEASAT  altimetry  (in  which 
they  allow  for  ample  redundancy  with  respect  to  the  ground- track  configuration), 
it  is  now  necessary  to  examine  all  the  altimeter  measurements  in  the  original 
observational  intervals  of  approximately  2*. 

The  most  important  part  of  the  design,  then,  resides  in  filling  in  the 
residuals  at  every  data  location.  This  task  can  be  accomplished  efficiently  upon 
undertaking  a  new  iteration  in  the  S.H.  adjustment  process,  but  this  time 
utilizing  all  the  altimeter  observations  and  limiting  the  iteration  merely  to  the 
computation  of  constant  terms  in  the  observation  equations.  These  terms  are  the 
desired  high-density  and  high-resolution  residuals  along  satellite  profiles. 

They  are  subsequently  judged  against  the  yardstick  provided  by  the  "truncation 
sigma"  of  the  underlying  S.H.  model,  i.e.,  the  square  root  of  the  variance  owing 
its  existence  to  the  truncation  of  theoretically  infinite  series.  As  an  example, 
a  sea  surface  dip  larger  than  two  or  three  truncation  sigmas  is  unlikely  to  re¬ 
flect  the  usual  geoidal  roughness,  and  may  instead  indicate  the  proximity  of  an 
important  trench.  Tests  in  areas  of  known  bathymetric  features,  especially  deep 
trenches,  constitute  a  useful  tool  in  a  "fine-tuning"  of  such  a  detection  pro¬ 
cedure. 

An  extensive  portion  of  the  present  research  has  been  devoted  to  attaining 
a  detailed  S.H.  representation  of  the  earth's  gravity  field  from  a  discrete  set 


of  data.  The  quantities  analyzed  most  thoroughly  in  this  context  are  the  geoid 
undulations,  directly  related  to  the  disturbing  potential  through  the  familiar 
Bruns  formula.  The  points  of  known  geoid  undulations  (and,  in  several  instances, 
gravity  anomalies)  are  considered  here  to  form  an  equilateral  grid.  The  main 
task  consists  in  finding  the  highest  degree  and  order  S.H.  model,  (N,N),  which 
can  still  accurately  represent  the  gravity  field  as  described  by  the  discrete 
data. 

This  task  is  addressed  through  computer  simulations  taking  advantage  of  the 
spherical  approximation.  As  a  first  step,  "errorless"  geoid  undulations  are 

generated  with  given  S.H.  coefficients  in  two  equilateral  grids,  referred  to  as 

0  0  0  0  0  0 

the  2  x2  and  the  4  x4  grid.  The  notation  6  x0  implies  that  most  of  the  grid 

lines  (except  near,  and  at,  the  poles)  give  rise  to  squares  with  0°  on  the  side. 
The  data  generating  process  hinges  on  an  algorithm  able  to  produce  accurate  values 
of  Legendre  polynomials  and  associated  Legendre  function  to  a  sufficiently  high 
degree  and  order.  One  such  algorithm,  extensively  tested  through  a  (160,160) 
truncated  model,  is  described  in  Appendix  1.  It  is  based  on  recurrent  relation¬ 
ships,  and  its  fundamental  feature  resides  in  utilizing  two  starting  values  for 
each  order  of  the  associated  Legendre  functions,  similar  to  the  strategy 
customarily  used  for  Legendre  polynomials.  A  great  increase  in  accuracy  compared 
to  an  approach  where  two  starting  values  are  used  to  generate  the  entire  set  of 
these  functions  is  well  worth  the  additional  effort. 

Linder  certain  conditions,  numerical  integration  involving  errorless  data 
leads  to  the  S.H.  coefficients  which  can  satisfactorily  reproduce  these  data 
at  the  grid  locations  as  well  as  at  other  points,  whether  in  their  original  form 
(here  geoid  undulations)  or  in  a  parent  form  (here  gravity  anomalies).  The 
theoretical  development  of  such  a  procedure,  together  with  all  the  pertinent 
computer  results  and  tests,  are  presented  in  Chapter  5.  The  key  element  in  this 


analysis  is  the  numerical  integration  algorithm  formulated  in  equations  (5.5a,b). 

There  is  no  need  to  recapitulate  much  of  the  material  of  Chapter  5  which  is 
straightforward  and  self-contained.  Its  most  important  outcome  is  represented 
by  Fig.  1,  depicting  the  errors  associated  with  the  numerical  integration  of 
geoid  undulations  in  the  above  20*2°  and  4°x4°  equilateral  grids.  An  error  measure 
is  provided  by  the  r.m.s.  of  the  differences  between  the  errorless  values  generated 
with  a  given  (n',n')  set  of  S.H.  coefficients  and  the  values  obtained  via  the 
integrated  coefficients  complete  to  within  the  same  set.  The  results  are  so  re¬ 
vealing  and  clearcut  that  the  curve  marked  "standard"  in  the  figure  becomes 
unnecessary.  This  curve  symbolizes  the  values  "t"  against  which  the  quality  of 
the  S.H.  representation  was  to  be  tested. 

Due  to  the  sharp  upturns  in  both  of  the  "2°x2°"  and  "4°x40"  curves,  occurring 
before  the  latter  intersect  the  "standard"  curve,  the  final  truncations  are  es¬ 
tablished  to  the  left  of  the  corresponding  intersections.  Thus,  the  final  outcome 
points  to  the  (90,90)  truncation  in  the  case  of  the  2°x2°  equilateral  grid,  and  to 
the  (45,45)  truncation  in  the  case  of  the  4°x4°  grid,  exactly  as  stipulated  by  the 
familiar  rule  of  thumb, 

N  =  180° /6°  . 

Clearly,  the  importance  of  this  rule  is  heightened  in  this  concrete,  tangible 
context. 

The  various  truncations  (n',n')  examined  in  Chapter  5  in  conjunction  with 
the  two  equilateral  grids  reveal  special  properties  worth  further  elaboration. 

As  Fig.  1  indicates,  the  " f* curve  is  effectively  a  straight  line  up  to  the 
truncation  (n' ,n‘ )=(90,90) ,  at  which  point  it  abruptly  changes  direction  and 
continues  again  as  a  straight  line,  but  at  a  much  greater  slope.  For  the  sake 


of  interest,  it  can  be  mentioned  that  a  curve  pertaining  to  gravity  anomalies, 
if  plotted,  would  exhibit  almost  identical  features.  The  "4°x4°"  <"  Jrve,  whose 
turning  point  is  (n‘ ,n‘ )=(45,45) ,  behaves  in  a  similar  fashion,  although  its 
lower  portion  departs  somewhat  from  a  straight  line  for  n‘  between  30  and  45. 

One  can  also  observe  from  the  figure  that  the  r.m.s.  of  the  discrepancies  in  the 
2°x2°  case  is  approximately  one-half  of  that  in  the  4° *4°  case,  at  least  up  to 
n'=30.  This  improvement  is  imputable  to  the  two-foid  increase  in  density,  along 
either  dimension,  of  the  former  grid  with  respect  to  the  latter. 

As  is  apparent  already  from  Fig.  1,  the  final  truncation  values  ( N ,N) =( 90 , 90) 
and  (N ,N )=( 45 ,45)  presented  above,  and  any  truncations  below  such  levels,  lead  to 
a  very  satisfactory  reproduction  of  the  errorless  geoid  undulations.  This  con¬ 
clusion  can  be  illustrated  quantitatively  upon  combining  (in  the  quadratic  sense) 
the  original  truncation  sigma  with  the  r.m.s.  of  the  discrepancies  due  to  the 
numerical  integration  algorithm,  and  obtaining  a  "new"  truncation  sigma  for  geoid 
undulations.  In  either  of  the  two  cases  under  consideration,  the  truncation 
sigma  is  seen  to  increase  by  a  mere  1%  --  an  encouraging  sign  indeed. 

The  outcome  of  Chapter  5  can  be  effectively  exploited  in  representing  colloca¬ 
tion  results  in  terms  of  the  S.H.  expansion  of  the  geopotential.  This  topic  is 
treated  in  Chapter  4  and  the  first  section  of  Chapter  6.  For  example,  a  2°*20 
equilateral  grid  of  geoid  undulations  predicted  through  the  modified  collocation 
with  noise  at  the  smoothing  level  (n1  ,n')  =  (90,90)  can  be  utilized  to  produce  an 
integrated  (90,90)  set  of  S.H.  coefficients  consistent  with  the  collocation 
results.  Since  the  emphasis  in  the  final  representation  is  on  the  oceanic  geoid, 
the  (missing)  prediction  values  over  land  can  be  filled  in  with  a  lesser  accuracy, 
upon  using  an  a  priori  set  of  S.H.  potential  coefficients  truncated  at  the  same 
level . 


Finally,  the  second  section  of  Chapter  6  addresses  the  problem  of  using 
a  S.H.  expansion  in  representing  the  desired  tidal  effects.  These  effects  are 
described  in  detail,  and  to  a  high  level  of  accuracy,  by  E.  w.  Schwiderski  in 
his  NSWC  Technical  Reports.  Each  tidal  constituent  considered  involves  two 
sets  of  S.H.  tidal  coefficients.  Some  of  the  principles  presented  in  Chapter 
5  are  utilized  here  as  well,  especially  the  numerical  integration  algorithm 
recapitulated  as  (6.3a-c).  This  algorithm  is  used  in  conjunction  with  an 
equilateral  grid  of  functions  called  F  and  G,  evaluated  over  water  with  the  aid 
of  the  above  NSWC  reports  and  over  land  via  visual  interpolation. 

If  this  grid  is  chosen  as  4° x4° ,  the  earlier  rule  of  thumb  allows  for  a 
reliable  reproduction  of  the  functions  F  and  G  through  a  (45,45)  set  of  integrated 
S.H.  tidal  coefficients.  A  lower  degree  and  order  set  can  also  be  used  in 
practice,  implying  further  smoothing  of  the  NSWC  information  (in  the  above  reports, 
the  tidal  information  is  available  in  a  l°xl°  geographical  grid).  Finally,  the 
ocean  tide  for  a  desired  constituent  “j"  can  be  expressed  as  in  (6.4),  and  can 
serve  in  the  tidal  adjustment  in  exactly  the  same  manner  as  its  former  counterpart 
did  in  CBlaha,  1982]  and  Chapter  3  of  CBlaha,  1984]. 


APPENDIX  1 


DEVELOPMENT  OF  A  RECURRENT  ALGORITHM 
FOR  LEGENDRE  POLYNOMIALS  AND  ASSOCIATED  LEGENDRE  FUNCTIONS 

This  Appendix  describes  in  detail  the  algorithm  developed  and  tested  during 

the  period  July  1  -  September  30,  1982,  and  presented  in  Status  Report  No.  3  for 

the  same  period.  A  similar  algorithm  was  summarized  in  September,  1982  by  Rapp 

C1982],  pages  9  and  10.  There  are  a  few  minor  differences  between  the  two 

approaches,  however.  For  example,  here  the  recurrent  relationships  are  derived 

in  terms  of  the  quantities  Ptm'1  rather  than  P  .  The  derivations  in  this  Appendix 

n  nm 

are  presented  against  the  AFGL  background  including  earlier  versions,  and  have 
thus  a  somewhat  tutorial  character. 

In  the  past,  geoid  undulations  and  other  geophysical  quantities  were  computed 
at  AFGL  on  a  global  scale  using  the  conventional  spherical-harmonic  (S.H.)  coef¬ 
ficients,  as  well  as  the  conventional  Legendre  polynomials  P  (sin<p)  and  the  con- 

n 

ventional  associated  Legendre  functions  Pnm(sin<j>),  where  <j>  denotes  the  geocentric 
latitude.  With  regard  to  the  Legendre  functions,  the  following  transformations 
apply: 

Pn(sin4>)  =  (2n  +  l)*5  Pn(sin<|>)  ,  (Al.la) 

Pnn/Sin<^  =  C2(2n  +  i) (n  -  m)  !/{n  +  m)  H^P^sin^)  ,  m  >  1  .  (Al.lb) 

The  products  between  the  S.H.  coefficients  and  the  corresponding  Legendre  functions 
are  the  sane  whether  the  conventional  (unbarred)  or  the  normalized  (barred) 
quantities  are  utilized. 

In  using  the  formula 

Pnm( Sin4>)  =  C0Sm<!’  Pn"':i(sin<l))  »  (A1.2) 
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where,  for  m  <  n  , 


PCm3(sin*)  =  (A5  (sin4>)/d(sin4>)m  (A1.3) 

n  n 

(for  m>n  this  expression  is  zero),  the  following  recurrent  relationships  have  been 
utilized  at  AFGL,  in  which  P  (sin4>),  P^^sin*)  are  written  simply  as  P  ,  PM  : 

Pn  =  (l/n)Csin^(2n- l)Pn  l- (n- l)Pn  23  ,  (A1.4) 

PCm]  =  PCm3  +  (2n- l)PCm:13  ,  2  <  m  <  n  -  2  ;  (A1.5) 

n  n-  2  n- 1  —  —  '  ' 

for  the  special  cases  m  =  1 ,  m  =  n-l,  m  =  n  the  formula  (A1.5)  becomes 


Cl3  . 
n 

(m=l) 

(A1.6a) 

Cn-1] 

=  (2n-  l)P^n:23  , 

(m=n  -  1) 

(Al. 6b) 

n 

n-l 

Cn3  = 

n 

(2n-  l)PCn"13  . 

n-l 

(m=n) 

(A1.6c) 

The  case  m  =  0  is  treated  in  (A1.4)  since  PnQ  =  P^  =  Pn  (see  equations  Al.2,3). 
The  two  values  needed  to  start  the  process  (A1.4)  are 

PQ=  1  ,  Px  =  sin*  . 

With  regard  to  the  subsequent  computation  of  pjjm3  ,  the  case  m 
according  to  (A1.6a).  The  starting  values  are 

(p£13  =  0)  ,  Pj13  =  1  .  (A1.8) 


(A1.7) 

= 1  is  treated  first 


The  first  of  these  values  is  trivial  (see  equation  A1.3)  and  is  not  stored;  it  is 
only  used  once  in  (A1.6a)  for  n = 2.  Next,  (A1.5)  is  utilized  for  each  n  (starting 
with  n  =  2)  and  for  each  m  (starting  with  m  =  2)  except  for  the  last  two  values  of 


m,  in  which  case  (A1.6b,c)  are  used  instead.  Thus  for  n  =  2,  only  (A1.6c)  is  used 
in  computing  P^23;  for  n  =  3,  (A1.6b)  gives  p|j23  and  (A1.6c)  gives  P^33;  but  for 
any  further  functions,  all  of  (A1.5)  (A1.6b,c)  are  needed. 

For  large  n  and  m,  the  values  of  P^  and,  especially,  P^m3  computed  as  above 
become  exceedingly  large.  Already  for  n  = 38  values  emerge  whose  order  of  magnitude 
is  lO70  and  higher,  leading  to  numerical  difficulties  because  of  the  computer 
exponent  limitations.  The  above  formulas  could  safely  be  used  to  within  a  (36,36) 
expansion.  Beyond  that,  the  normalized  functions  should  be  used  instead.  New 
recurrent  relations  can  be  obtained  from  the  previous  formulas  through  a  simple 
substitution  for  P  ,  P1"™3  based  on  (Al.la.b),  namely 

n  n 

P  =P/(2n+l)%,  (A1.9a) 

n  n  ' 

P^  =  {(n  +  m)!/C2(2n+1)(n-m)!3}Vnm  ,  m  >  1  .  (A1.9b) 

Since,  in  analogy  to  (A1.2),  it  also  holds  true  that 

P  =  cos%  PCm3  ,  (A1.10) 

nm  n 

the  formulas  (Al.lb)  and  (A1.9b)  relate  P^m3  and  P^m3  as  well.  This  will  be 
understood  throughout  without  any  further  statements  to  that  effect. 

Equation  (A1.4)  can  now  be  rewritten  with  (A1.9a)  taken  into  account,  giving 
P  =  (2n  + l)l5(l/n){sin^(2n- D^P  -  c(n  -  l)/(2n  -  3)SDP  ,}  .  (Al.ll) 

n  n-1  n-Z 

If  (A1.5)  -  (A1.6c)  are  similarly  rewritten  with  (A1.9a,b)  taken  into  account,  it 
follows  that 

P^m3  =  {(2n  +  l)/[(n+m- l)(n +  m)3  }S{C(n  -  m  -  1 )  (n  -  m)/(2n  -  3)]V;™3 

+  (2n  -  l)VCm'13}  ,  2  <  m  <  n  -  2  ;  (A1.12) 

n-l  —  — 
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P1-13  =  {(2n  +  l)/Cn(n  +  l)D}1,{C{n-2)(n-l)/(2n-3)D1,PI:i3  +  C2(2n-l)D‘IP  ,}  , 

n  n-  2.  n- 1 

(m  = 1)  (Al. 13a) 

PCn"13  =  C(2n  +  l)/(2n- 2)3lsPCn723  ,  (m  =  n-l)  (Al.  13b) 

n  n-1 

P*”3  =  c(2n  +  l)/2n31*P^n713  .  (m  =  n)  (Al .  13c) 

n  n-1 

Due  to  the  presence  of  2  in  the  factor  of  P  in  (A1.13a),  this  formula  is  not  a 
special  case  of  (A1.12),  in  the  sense  that  it  cannot  be  obtained  upon  replacing 
m  in  (A1.12)  by  1  (this  could  have  been  done  when  passing  from  A1.5  to  A1.6a). 

In  agreement  with  (A1.7)  and  (Al.la),  the  starting  values  for  (Al.ll)  are 

P Q  =  1  ,  P1  =>/3  sin*  .  (Al. 14) 

Subsequently,  in  using  ( Al. lb)  the  values  of  P^m3  are  computed  in  analogy  to 
(A1.8)  and  the  text  that  followed;  the  starting  values  are 

(Pq13  =  0)  ,  P*13  =  >/J  .  (Al. 15) 

Large  numerical  values  are  no  longer  a  hindrance  in  proceeding  in  this  fashion. 
However,  for  large  n  and  m  --  around  70  or  80  and  higher  —  the  computational 
accuracy  in  F^m3  quickly  deteriorates.  This  is  due  to  the  propagation  of  round¬ 
off  errors  since  a  great  many  values  (on  the  order  of  3,000  and  higher)  are 
generated  from  merely  two  starting  values  in  (A1.15).  Therefore,  the  present 
set  of  formulas  could  safely  be  used  only  to  about  a  (70,70)  expansion.  Although 
this  is  a  significant  improvement  compared  to  the  previous  situation,  it  does  not 
fulfill  the  anticipated  needs.  Consequently,  a  new  set  of  recurrent  relations 
will  be  derived  which  can  be  used  for  any  expansion  in  practice,  certainly  thrjugh 
the  degree  and  order  (180,180)  and  possibly  much  higher. 


The  basic  idea  behind  the  new  approach  is  to  use  two  starting  values  for 
each  order  (m)  of  the  associated  Legendre  functions,  similar  to  the  recurrent 
strategy  for  Legendre  polynomials.  Thus  in  the  case  of  the  truncation  N  =  180, 
at  most  178  values  will  be  generated  from  two  starting  values.  In  particular, 
since  one  will  proceed  in  the  sequence  PCm3,  P^m3 . PM,  the  number  of 

m  m+1  n 

generated  values  varies  from  178  for  m  =  l  (P^13  and  P^13  are  the  starting  values 
while  P^1],  P*;13,  ...  ,  are  generated)  to  1  for  m  =  178  (P^I83  and  P^I8] 

j  4  180  i/o  l/y 

_ r  1 7  qt 

are  the  starting  values  while  only  Pr  is  generated);  no  values  are  generated 

180 

for  m  =  179  (both  P*^93  and  P^I93  are  "starting"  values)  and  for  m  =  180 

179  180 

(P^803  is  a  "starting"  value).  Although  there  are  essentially  N-times  as  many 
starting  values  in  this  approach  compared  to  that  of  the  previous  paragraphs, 
the  great  increase  in  accuracy  is  well  worth  the  effort. 

This  development  starts  with  the  conventional  Legendre  functions.  First, 
the  relation  (A. 14)  for  Legendre  polynomials  is  recapitulated  as 


nP  =  s i net  ( 2n  -  1 ) P  ,  -  (n-l)P  _  ,  n  >  2  . 

n  n-1  n-2  — 


(A1.16) 


This  formula  can  be  found  in  standard  textbooks  such  as  CSpiegel ,  196811  where  it 
is  essentially  equation  (25.20).  The  starting  values  are  PQ  and  P1  as  they 
appeared  in  (A1.7).  Next,  the  basic  relation  for  P^m3  is  presented  as 

(n-rn)PCm3  =  sin$  (2n- l)PCm3  -  (n  +  m-  l)PCln3  ,  n>m+2,  (A1.17) 

n  n-1  n- 2  — 


which  is  readily  obtained  from  equation  (26.12)  in  the  same  textbook.  This  formula 

is  structured  in  complete  analogy  to  (A1.16),  which  can  be  obtained  from  it  merely 

by  setting  m  =  0.  The  starting  values,  P M  and  PCm3  ,  can  be  evaluated  with  the 

m  m+1 

aid  of  the  standard  formula  appearing,  e.g.,  as  (1-62)  in  CHeiskanen  and  Moritz, 
i9673,  which  leads  to 
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(Al. 18a) 


ptmj  =  (1/2m)(2m)  !/m!  =  1x3  ...  *  (2m- 1)  , 
m 

=  1x3  ...  x  (2m  +  l)sin4>  .  (A1.18b) 

m+1 

Here  again,  upon  setting  m  =  0  one  recovers  the  starting  values  (A1.7). 

In  utilizing  (A1.9a)  in  (A1.16)  it  follows  that 

Pn  =  (2n  +  D^d/n)  (sin4»  (2n-l)l*Pn_1  -  C(n  -  l)/(2n  -  3)^3Prv_2)  ,  (A1.19) 

which  is  (Al. 11)  recapitulated  here  to  show  the  parallel  between  the  new  recurrent 
relations  for  and  for  P^m]  below.  In  considering  (Al.la),  the  overbarred 
starting  values  are 


with  PQ  and  P1  from  (A1.7)  this  yields 

P 0  =  1  ,  (A1.20a) 

?1  =  sin#  ,  (Al .  20b) 

which  were  already  seen  in  (A1.14).  Equation  (Al .19)  with  the  starting  values 
(A1.20a,b)  is  the  final  recurrent  relation  for  the  normalized  Legendre  poly¬ 
nomials. 

In  analogy  to  the  above,  upon  utilizing  (A1.9b)  in  (A1.17),  after  a  few 
arithmetic  operations  one  obtains 

P^1"3  =  {(2n  +  l)/c(n  +  m)(n  -  m)]  }!3{sin4>  (2n  -  l)^1"3  -C(n  +  m-  l)(n-m-l) 

n  n- 1 

/(2n- 3)Dl5PDn:'}  ,  n  >  m  +  2  .  (Al. 21) 

n-2 
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a'*  A 


>  >  »v.v: 


It  is  again  noticed  that  upon  setting  m  =  0  equation  (A1.19)  is  recovered.  In 
considering  (Al.lb),  the  overbarred  starting  values  are 


pCm:  =  £2(2m  +  l)/(2m)  ,  m  >  1  , 

m  m  — 

P1-1"?  =  C2(2m  +  3)/(2m  + 1)  ,  m  >  1  ; 

m+i  m+1  — 

with  PCn,:i  and  PCm?  from  (A1.18a,b)  this  yields 
m  m+1 

P1"13  =  C2(2m  +  l)3*5{(l/2)  (3/4) . . .  C(2m  -  l)/2m3  J*3  ,  m  >  1  ,  (A1.22a) 

m  — 

pk"?  =  C2(2m  +  3)]ls{{3/2) (5/4) . . . c(2m  +  l)/2m] l^sin^  ,  m  >  1  .  (Al- 22b) 

m+1  — 


It  is  now  apparent  that  (A1.20a,b)  do  not  follow  from  (A1.22a,b)  upon  setting 
m  =  0,  due  to  the  factor  2  within  the  first  brackets  of  the  latter.  This  is 
imputable  to  the  factor  2  appearing  in  equations  (Al.lb)  and  (A1.9b)  which  thus 
would  not  yield  (Al.la)  and  (A1.9a)  if  m  were  set  to  zero;  this  fact  is  then 
reflected  in  the  equations  that  followed  (A1.21)  as  compared  to  the  equations  that 
followed  (A1.19).  On  the  other  hand,  (A1.19)  followed  from  (A1.21)  for  m  =  0, 
just  as  (A1.16)  followed  from  (Al. 17)  for  such  m,  because  when  the  transformation 
(A1.9b)  was  applied  to  (A1.17)  in  order  to  obtain  (A1.21),  the  above  factor  2 
cancelled  out. 


Equation  (A1.21)  with  the  starting  values  (A1.22a,b)  is  the  final  recurrent 
The  final  recurrent  relation  for  the  normalized  associated 
Legendre  functions  then  follows  immediately  upon  considering  (A1.10).  In  parti- 


relation  for  P^ 

n 


cular,  since  no  order  other  than  m  appears  in  (A1.21)  this  equation  holds  true 
equally  well  with  P  ,  P  .  and  P  ,  replacing  F?"'3,  P*^  and  P^  , 
respectively.  Since  the  starting  values  are  given  by  the  relatively  very  simple 
formulas  (A1.22a,b),  where  the  number  of  factors  within  (1  is  m,  their  computa¬ 
tion  does  not  represent  a  heavy  burden  by  any  standards.  In  fact,  without  the 


APPENDIX  2 


CONSIDERATION  OF  SPECIAL  "TIDAL  POINT  MASSES" 

FOR  MODELING  OF  TIDAL  CONSTITUENTS 

In  CBlaha,  1984],  Chapter  3,  the  tidal  model  used  in  the  adjustment  of 
satellite  altimetry  is  given  as 

lv  =  0.9333^  +  0.612x(equil ibrium  tide)..  ,  (A2.1) 

where 

h  =  surface  tide  for  the  constituent  "j"  , 

3 

=  ocean  tide  for  the  constituent  "j"  . 

The  formula  (A2.1),  appearing  in  the  above  reference  as  (3.7a,b),  can  be  used  in 
conjunction  with  all  the  tidal  constituents,  'c.ig-period  as  well  as  diurnal  and 
semidiurnal . 

Since  the  computation  of  the  equilibrium  tide  is  simple  and  straightforward, 
one's  attention  is  directed  to  the  evaluation  of  the  ocean  tide  in  equation  (A2.1) 
According  to  pages  48  and  49  in  CBlaha,  1982],  it  follows  that 

^  =  a  ((jl,x)coscoj  -  ^U.x)]  , 

where 

a  =  Greenwich  argument, 

A  ($»x)  =  amplitude, 

(<l>,x)  =  Greenwich  phase  angle, 

all  pertaining  to  toe  constituent  "j". 


The  above  equation  can  be  written  as 


with 


(A2.; 


£  =  COSa  .A .  ($>,a)cosi|j  .  ( 4> » A )  ,  (A2.I 

3  J  J  J 

£  =  sina.A.(i}i,A)sin^.((j),A)  .  (A2.I 

3  3  3 

The  functions  F  and  6  below  can  be  expressed  in  terms  of  a  spherical-harmonic 
series  as  was  done  in  the  above  reference,  namely 


F(*,a)  5  AjU,A)cos<|/jU,A)  = 
G(<M)  =  Aj  (4>,A)sin^j(4>.A)  = 


(A2. 

oo  n 

z  z  (a.  cosmA  +  b.  sin  mA)P  (sin<j>)  , 

n=0  m=0  -'nm  -)nm  nm 

(A2. 

oo  n 

z  z  (c.  cos  mA  +  d.  sin  mA)P  (si n<j>)  , 

*  _  3  3  ran 

n=0  m=0  ran  nm 


where  P  ( sincfr)  are  the  associated  Legendre  functions  in  the  argument  s i n<j>  , 

4>  and  a  being  the  geocentric  latitude  and  longitude,  respectively. 

If  the  functions  F  and  G  were  given  all  over  the  globe,  one  could  compute 
the  spherical-harmonic  (S.H.)  tidal  coefficients  a,  b,  c,  d  (with  appropriate 
indices)  as  follows: 


a . 

-'ran 

b. 

^  nm 


C(2n+l)/2TT]C(n-m)  !/(n+m)  !D//cA.(^,a)cos^ . ($, x) DP  (sin$) 
—  0  3  3  ran 


COS  mA 
sin  mx 


c . 
-'ran 

A 


cos  mx 


C(2n+l)/2Tr3t(n-m)  !/(n+m)  !3//tA .  ( ^ ,  x)  si  nip  .(^,a)DP  (sin^) 


where  the  symbol  2n  .mplies  2tt  for  m>0  and  4n  for  m=0;  this  double  role  would 
have  been  avoided  by  adopting  the  normalized  S.H.  coefficients  and  Legendre 


functions  instead  of  the  conventional  ones  as  used  above.  The  symbol 
indicates  a  point  within  do,  a  representing  the  surface  of  the  earth. 

With  the  expansions  (A2.3a,b),  equations  (A2.2a-c)  yield 

oo  n 


£  . 


3 


COSa  . 
3 


i  z  (a.  cos  mx  +  b.  sin  mx)P  (sin<j>) 

n=0  m=0  ^nm  ^nm  riri 


oo  n 


+  sina. 
3 


z  z  (c.  cos  mx  +  d .  sin  mx)P  (sin<|>)  . 

n=0  m=0  ^nm  ^nm  lm 


II  ->  II 


(A2.4) 


The  parts  following  cosou  and  sina^  in  (A2.4)  can  be  compared  to  the  S.H. 
expansion  of  geoid  undulations  (N).  However,  it  is  assumed  that  a  first-phase 
adjustment  has  been  carried  out  so  that  the  second-phase  adjustment  has  as  a 
"normal  field"  the  gravity  field  given  by  the  (n,n)  S.H.  expansion.  With  respect 
to  such  a  field,  the  second-phase  geoid  undulations  resolved  through  the  degree 
n'  are  represented  by 

n’  k 

N'  =  R  z  Z  {C,  cos  mx  +  S,  sin  mx)P,  (sinij>)  ,  (A2.5) 

,  ,  „  Km  Km  Km 


where  R  is  the  earth's  mean  radius  and  <t>,x  are  the  geocentric  coordinates  of  the 
point  of  evaluation.  The  value  of  N'  in  (A2.5)  thus  corresponds  to  the  bandwidth 
n+1  through  n'.  Since  n  is  usually  14  or  higher,  the  spherical  approximation 
implied  by  (A2.5)  is  inconsequential. 

The  similarity  between  (A2.4)  and  (A2.5)  is  actually  closer  than  what  appears 
from  the  equations  in  their  current  form,  especially  from  the  form  of  (A2.4).  The 
first-phase  adjustment  provides  not  only  for  the  solution  of  the  gravity  field 
within  the  (n,n)  S.H.  expansion,  but  also  for  the  solution  of  selected  tidal 


parameters  to  within  some  expansion  of  the  S.H.  tidal  coefficients.  The  tidal 
parameters  consist  of  the  change  in  the  relative  amplitude  and  the  change  in 
the  phase  angle  per  tidal  constituent;  the  degree  and  order  of  the  underlying 
S.H.  expansion  of  the  tidal  effects  in  the  first  phase  is  considered  to  be  also 
(n,n),  although  this  is  not  a  requirement  (the  advantage  of  this  choice  stems 
from  the  Legendre  functions  having  been  already  computed  to  within  n,n  degree 
and  order  due  to  the  geoidal  evaluations,  and  the  same  is  true  with  respect  to 
the  multiple  angle  trigonometric  functions).  Accordingly,  if  one  wished  for 
the  solution  of  the  ocean  tide  within  the  bandwidth  n+1  through  n'  as  a  follow¬ 
up  to  the  first  phase,  the  formula  (A2.4)  for  a  given  tidal  constituent  (the 
explicit  subscript  j  is  dropped)  would  become: 

n'  k 

£'  =  cosa  e  e  (a,  cos  mx  +  b,  sin  mx)P,  ( s i n<j> ) 

k=n+l  m=0  km 

(A2.6) 

n'  k 

+  sina  E  E  (c.  cos  mA  +  d,  sin  mx)P,  (sin<j>)  . 

.  -  ~  Km  Jan  km 

k=n+l  m=0 

The  complete  analogy  between  the  formula  (A2.5)  and  either  part  of  the  formula 
(A2.6)  is  apparent. 

The  above  analogy  suggests  that  the  second-phase  treatment  of  tidal  param¬ 
eters  could  follow  the  lines  similar  to  the  second-phase  treatment  of  gravity 
field  parameters.  The  latter  have  been  represented  by  point-mass  (P.M. )  magni¬ 
tudes  allowing  the  local  effects  to  manifest  themselves  quite  independently 
of  other  such  effects  some  distance  apart,  due  in  part  to  the  cut-off  distance 
(spherical  cap)  beyond  which  the  contribution  of  observations  to  a  given  P.M. 
parameter  has  been  ignored  by  the  algorithm.  Accordingly,  the  tidal  effects 
could  also  have  a  greater  freedom  to  manifest  themselves  locally  if  parallel 


formulas  could  be  developed  in  terms  of  some  "tidal  point  masses".  The  second- 
phase  adjustment  model  in  terms  of  the  “regular"  point  masses  reads 

=  (1/G)  E  (l/*kj)(kM) .  ,  (A2.7) 

where  G  is  the  average  value  of  gravity  on  the  earth's  surface,  "k"  is  the 
(observation)  point  where  the  geoid  undulation  is  available,  "j"  is  the  location 
of  the  point  mass  whose  magnitude  is  (kM)^,  and  j,  is  the  (chord)  distance 
between  the  points  k  and  j.  Since  (A2.7)  is  used  to  approximate  the  model  (A2.5), 
the  analogy  between  (A2.5)  and  (A2.6)  allows  the  latter  to  be  similarly 
approximated  by 

^  =  cosa  i  (1/^.)^  +  sina  E  (l/fck.)t2  .  (A2.8) 

j  3  j  j  3  j 

Here  a  remains  the  Greenwich  argument  and  ^  ,  t2  represent  the  "tidal  P.M. 

j  j 

parameters" ,  all  pertaining  to  a  given  tidal  constituent. 

Since  the  geoid  undulations  and  the  effects  of  tidal  constituents  are 
mutually  independent,  so  are  the  P.M.  parameters  and  the  tidal  P.M.  parameters 
for  different  constituents.  The  latter  can  thus  be  added  into  the  system  of 
normal  equations  following  the  former  without  affecting  the  bandwidth  associated 
with  a  "regular"  P.M.  adjustment.  However,  the  length  of  the  band  is  greatly 
increased.  If  the  tidal  point  masses  should  have  the  same  locations  as  the 
"regular"  point  masses,  the  number  of  parameters  would  increase  by  200%  per 
constituent  as  is  clear  upon  comparing  the  model  equations  (A2.7)  and  (A2.8). 

Where  only  the  constituents  M2,  N2  and  0  are  considered,  the  increase  would 
be  600%,  so  that  the  total  number  of  parameters  would  be  7 *(" regular"  number  of 
P.M.  parameters).  Such  an  increase  would  be  clearly  prohibitive. 


However,  the  tidal  constituents  have  a  much  smaller  effect  on  altimeter 
residuals  than  the  geoid  undulations  (one  to  two  orders  of  magnitude). 

Accordingly,  a  low  density  grid  of  tidal  point  masses  would  be  sufficient  for  a 
simultaneous  adjustment  of  geoidal  and  tidal  effects.  For  example,  one  can  con¬ 
sider  a  2°-geoidal  resolution  corresponding  approximately  to  a  (90,90)  S.H. 
expansion  of  the  earth's  gravity  field  contrasted  to  a  6°-tidal  resolution  cor¬ 
responding  approximately  to  a  (30,30)  S.H.  expansion;  the  number  of  S.H.  coef¬ 
ficients  in  the  latter  case  is  about  nine- fold  smaller  than  in  the  former  case. 

The  tidal  resolution  could  have  been  lowered  even  further  with  respect  to  the 
geoidal  resolution.  Be  that  as  it  may,  the  number  of  tidal  point  masses  is  about 
1/9  of  the  "regular"  point  masses;  when  multiplying  this  number  by  6  (two  param¬ 
eters  per  constituent,  three  constituents),  the  increase  in  the  number  of  param¬ 
eters  is  about  67%,  for  a  total  of  1.67x("regular"  number  of  parameters). 

The  dimensions  in  the  P.M.  algorithm  have  been  recently  increased  to  a  maximum 
of  2000  P.M.  parameters  by  virtue  of  the  banded  structure  of  normal  equations. 
However,  almost  that  many  are  needed  for  the  P.M.  adjustment  in  the  Indian  Ocean 
alone,  i.e.,  about  1700  P.M.  parameters  including  those  that  have  to  be  used  twice 
due  to  overlapping  adjustment  strips.  Multiplying  this  number  by  1.67  would 
result  in  a  total  number  of  parameters  reaching  about  2800.  This  number  of 
parameters  is  still  too  great  of  a  burden  regarding  the  computer  capacity. 

Next  consider  a  8°-tidal  resolution  corresponding  approximately  to  a  (22,22) 
S.H.  expansion.  The  tidal  point  masses  would  now  be  distributed  in  an  S'xS" 
equilateral  grid  and  their  number  would  be  1/16  of  the  "regular"  P.M.  parameters. 
The  multiplication  by  6  results  in  an  increase  of  37.5%  in  the  number  of  param¬ 
eters  for  a  total  of  1.375x("regular"  number  of  parameters);  in  the  above  case 
of  the  Indian  Ocean  this  implies  some  2300  parameters.  From  the  computer  core 
standpoint  an  adjustment  could  become  feasible.  The  algorithm  would  be  similar 


to  the  "regular"  P.M.  algorithm  except  for  some  self-evident  modifications  (see 
cosa  and  sinu  in  equation  A2.8,  the  introduction  of  two  parameters  per  location 
when  treating  a  given  constituent,  etc.)*  The  depth  of  the  tidal  point  masses 
could  again  correspond  to  0.8  of  their  horizontal  separation  (s);  this  would 
lead,  similar  to  previous  studies,  to  considering  a  spherical  cap  of  the  radius 
1.5s  when  determining  which  observations  should  be  ignored  in  conjunction  with 
a  given  parameter,  and  would  preserve  the  original  bandwidth  essentially  unchanged. 

However,  from  the  practical  point  of  view  the  above  case  would  not  have  the 
desirable  effect.  The  first-phase  adjustment  at  AFGL  is  usually  made  in  terms 
of  a  (14,14)  S.H.  expansion  of  the  earth's  potential,  and  a  higher  degree  and 
order  expansion,  within  perhaps  a  (20,20)  model,  is  easily  feasible.  Due  to  the 
advantages  already  mentioned,  the  S.H.  tidal  expansion  should  be  made  to  about 
the  same  degree  and  order.  Thus,  in  view  of  a  (20,20)  S.H.  expansion  of  tidal 
effects  being  essentially  a  part  of  the  first-phase  algorithm  already,  the  above 
(22,22)  expansion  is  of  marginal  usefulness.  In  other  words,  the  computer  core 
capacity  would  have  to  nearly  double  before  tidal  P.M.  parameters  in  the  second 
adjustment  offered  a  tangible  advantage. 
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APPENDIX  3 

UTILIZATION  OF  SEASAT  ALTIMETER  RESIDUALS 
IN  DETECTING  BATHYMETRIC  ANOMALIES 


As  has  been  reported  on  numerous  occasions,  pronounced  sea-bottom  features 
are  accompanied,  at  the  sea-level  surface,  by  localized  roughness  in  the  oceanic 
geoid.  This  fact  is  best  illustrated  by  the  presence  of  large  negative  alti¬ 
meter  residuals  in  the  area  where  a  satellite  arc  crosses  an  important  trench, 
such  as  the  Puerto  Rico  trench,  the  New  Hebrides  trench,  the  Aleutian  trench 
and  others.  Such  residuals  were  presented,  for  example,  in  the  paper  "SEASAT 
Altimeter  Reductions  for  Detailed  Determinations  of  the  Oceanic  Geoid"  by 
G.  Hadgigeorge,  G.  Blaha  and  T.  P.  Rooney  at  the  1980  AGU  Fall  Meeting. 

A  situation  of  this  kind  is  schematically  depicted  in  Fig.  A,  where  the 
(14,14)  geoid  corresponds  to  a  global  spherical-harmonic  (S.H.)  adjustment  of 
satellite  altimetry.  The  (90,90)  geoid,  based  on  the  residuals  from  the  global 
adjustment,  represents  relatively  small  features,  down  to  those  whose  half¬ 
wavelength  is  approximately  2°.  This  fairly  detailed  resolution  can  be  achieved 
through  the  use  of  point  masses  (distributed  in  a  2°x2°  equilateral  grid), 
collocation  predictions  (in  a  2°x2°  equilateral  grid)  and  other  means.  The  two 
methods  just  mentioned  are  described  in  the  AFGL  report  CBlaha,  1984];  the 
first  method  alone  has  also  been  the  subject  of  several  earlier  AFGL  reports 
and  papers.  Clearly,  the  above  (14,14)  and  (90,90)  S.H.  expansions  could  be 
replaced  by  other  suitable  degree  and  order  expansions  describing  the  earth's 
gravity  field  to  within  the  desired  resolution. 

The  index  "o"  in  Fig.  A  indicates  "observed"  quantities,  which  are  the 
originally  measured  values  improved  through  the  orbital  and  tidal  adjustment 
(but  not  the  geoidal  adjustment).  The  orbital  adjustment  is  implied  in  the 


=  adjusted 

(14,14)  geoid 

=  adjusted 
(90,90)  geoid 

=  sea  surface 
(tidal  effects 
considered 
removed) 


H.»0«  -  "observed"  altimetry 

Ha  =  adjusted  (14,14)  altimetry 
N„0„  =  "observed"  geoid  undulation 
Na  =  adjusted  (14,14)  undulation 


Fig.  A 


Schematic  representation  of  satellite  altimetry  over  a  trench  area 
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figure  by  the  altimeter  measurement  associated  with  the  adjusted  --  as 
opposed  to  the  initially  given  --  satellite  arc.  And  the  tidal  adjustment  is 
implied  by  considering  the  sea  surface  as  free  from  tidal  variations  which  are 
expressed  by  separate  parameters.  The  superscript  "a"  indicates  adjusted 
quantities,  in  particular,  quantities  linked  to  the  (14,14)  geoid  as  obtained 
in  the  global  S.H.  adjustment.  Figure  2  of  the  AFGL  report  CBlaha,  1983] 
illustrates  the  "observed"  and  adjusted  quantities  in  more  detail.  In  the  present 
context  one  need  not  be  preoccupied  with  the  improvement  due  to  the  tidal  adjust¬ 
ment  since  it  would  amount,  in  open  ocean,  to  merely  a  few  decimeters.  By 
comparison,  the  altimeter  residuals  over  important  trench  areas  are  expected  to 
be  greater  by  one  or  two  orders  of  magnitude. 

Altimeter  residuals  will  be  denoted  here  by  the  symbol  v  .  These  quantities 
are  treated  as  measurements  in  the  subsequent  determinations  of  the  gravity  field 
detail  to  be  superimposed  on  the  (14,14)  S.H.  representation,  upon  using  the 
point-mass  or  collocation  techniques  as  mentioned  earlier.  Denoting  the  geoidal 
residuals  by  the  symbol  v  .  one  has 

N 


•  Ha  -  H„o„  , 

(A3. la) 

=  Na  -  N„o„  , 

(A3. lb) 

(A3. lc) 

As  depicted  in  Fig.  A,  a  trench  area  induces  large  negative  altimeter  residuals, 
vN  «  0  .  (A3. 2) 

These  values  are  to  be  compared  with  the  "one  sigma"  such  as 


which  characterizes  the  geoidal  detail  neglected  due  to  the  (14,14)  truncation 
of  an  infinite  S.H.  series.  The  description  and  estimation  of  this  quantity  can 
be  found  in  Section  4.3  of  CBlaha,  19843. 

The  three  trench  areas  mentioned  at  the  outset  are  100-200  km  wide  and  are 
associated  with  large  negative  altimeter  residuals  reaching  typically  between 
-10  m  and  -15  m  in  the  vicinity  of  the  trench  axes.  One  could  detect  sea  surface 
dips  on  this  scale  from  a  detailed  geoidal  map  representing  a  Irresolution  or 
better  (1°  corresponds  to  111  km  on  the  earth's  surface).  In  order  to  detect 
dip  areas  100  km  wide,  a  %°-resolution  would  be  desirable.  But  such  a  procedure 
has  a  significant  drawback  imputable  to  the  spatial  distribution  of  SEASAT 
altimeter  data.  The  ascending  and  descending  passes  form  approximately  a  1°  *1° 
equilateral  grid  near  the  equator,  providing  for  a  strong  2°-resolution  and  a 
marginal  l°-resolution,  but  making  any  finer  resolution  impossible  to  achieve. 

In  addition  to  its  marginal  geoidal  representation,  the  l°-resolution  would  still 
be  too  coarse  for  detecting  relatively  small  but  important  ocean  trench  areas. 

The  notion  of  coarseness  in  this  context  is  best  illustrated  by  considering 
a  2°-resolution  as  obtained  through  a  point-mass  (P.M. )  adjustment  of  the  oceanic 
geoid.  The  observational  density  on  SEASAT  arcs  is  just  under  2'  in  angular 
measure.  Clearly,  utilizing  all  of  the  altimeter  data  in  the  S.H.  adjustment 
and,  subsequently,  in  the  P.M.  adjustment  would  be  computationally  prohibitive. 
Furthermore,  little  would  be  gained  from  a  configuration  where  the  observational 
density  along  one  dimension  is  over  30  times  higher  than  the  density  along  the 
other  dimension.  A  practical  conclusion  has  been  reached  to  adopt  every  16th 
point  on  each  arc  as  selected  input  data  in  both  adjustments.  In  this  way,  the 
separation  between  measurements  along  tracks  becomes  h° ,  while  the  separation 
across  tracks  is  fixed  at  1°.  This  allows  for  a  sufficient  amount  of  spatial 
redundancy  in  the  P.M.  adjustment,  especially  in  view  of  a  2° -resolution 
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(a  1° -resolution  with  the  same  data  distribution  would  imply  50%  redundancy 
along  tracks  and  no  redundancy  across  tracks). 

The  smoothing  effect  of  a  2°-resolution  could  totally  or  partially 
obliterate  the  representation  of  a  sea  surface  dip  in  the  proximity  of  a  trench. 
This  has  been  noticed,  for  example,  in  the  paper  mentioned  at  the  outset.  The 
adjusted  geoid  in  an  extended  Puerto  Rico  trench  area  covered  with  a  2°*2° 
equilateral  grid  of  point  masses  did  not  exhibit  any  marked  dip  until  it  was 
saturated  with  additional  point  masses  forming  a  localized  l°xl°  equilateral 
grid.  It  is  noteworthy  that  the  geoidal  dip  associated  with  the  Puerto  Rico 
trench  is  pronounced  (on  the  order  of  15  m),  that  the  trench  area  is  relatively 
wide  (approximately  200  km)  and,  especially,  that  the  trench  location  is  known 
beforehand.  A  less  important  trench  might  have  gone  unnoticed  altogether. 

The  example  of  the  Puerto  Rico  trench  is  overly  optimistic,  in  that  not 
only  is  the  trench  area  wide  and  the  location  known,  but  also  that  the  observa¬ 
tional  interval  along  tracks  is  This  presents  little  computational  burden  by 
virtue  of  the  regional  extent  of  the  adjustment,  but  the  same  cannot  be  said 
with  regard  to  an  adjustment  of  an  entire  ocean  basin  or  the  entire  global  ocean. 
One  might  hesitate  to  engage  in  the  exercise  of  a  global  P.M.  adjustment  with 
quadruple  the  number  of  parameters  (in  view  of  the  1° -resolution)  and  double  the 
number  of  data  points  (in  view  of  data  intervals  along  tracks)  merely  for  the 
purpose  of  detecting  several  more,  but  by  no  means  all,  important  trench  areas. 

However,  even  if  one  should  choose  to  carry  out  such  an  exercise  one  would 
be  well  advised  to  consider  the  adjusted  1° -geoid  together  with  the  resulting 
(second-phase)  residuals,  especially  if  the  detection  of  sea  surface  anomalies 
should  not  be  limited  to  features  whose  width  surpasses  100  km.  But  the  adjusted 
P.M.  geoid  plus  the  corresponding  residual  equal  the  original  altimeter  residual 
from  the  global  S.H.  adjustment.  This  is  apparent  from  Fig.  A,  where  the 


distance  between  the  solid  and  the  dotted  curve  is  obtained  by  adding  algebraical¬ 
ly  the  two  segments  created  by  the  dashed  curve.  Thus,  the  second-phase 
techniques,  such  as  the  P.M.  adjustment,  collocation  predictions,  etc.,  can  be 
completely  by-passed  in  the  present  task. 

The  use  of  the  first-phase  technique,  i.e.,  the  global  S.H.  altimeter  ad¬ 
justment,  has  a  further  significant  advantage.  In  particular,  a  procedure  has 
been  designed  that  takes  into  account  every  observation  point  (in  approximately 
2'  intervals  along  tracks),  not  only  points  in  *s°  or,  perhaps,  %°  intervals. 
Although  the  S.H.  adjustment  itself  may  utilize  only  data  in  intervals,  the 
remaining  residuals  can  subsequently  be  filled  in  at  every  data  location.  These 
high-density  and  high-resolution  residuals  can  be  examined  for  geoidal  roughness 
along  the  profiles,  which  has  a  crucial  role  to  play  in  detecting  bathymetric  and 
possibly  other  anomalies,  such  as  geomagnetic  anomalies. 

The  most  promising  in  this  respect  is  the  detection  of  ocean  trenches  whose 
presence  is  signalled  by  an  abrupt  dip  in  the  sea  surface.  The  importance  of  any 
such  dip  is  judged  upon  comparing  the  (first-phase)  altimeter  residuals  with  the 
value  a  linked  to  the  truncation  of  the  S.H.  model.  To  be  of  any  practical 
advantage  this  value  should  be  fairly  low  and,  accordingly,  the  S.H.  model  should 
be  of  a  reasonably  high  degree  and  order.  Equation  (A3. 3)  indicates  that  the 
(14,14)  model  is  indeed  satisfactory,  considering  that  altimeter  residuals  at  the 
dip's  location  reach  three  to  five  times  the  magnitude  of  a. 

The  satellite  arcs  entering  the  S.H.  adjustment  are  limited  in  length  to 
about  25°  in  angular  measure  or  7  minutes  in  duration  as  explained  in  Chapter  2 
of  CBlaha,  19813.  With  the  selected  data  points,  a  25°-arc  contains  54  obser¬ 
vations  in  ^-intervals,  as  opposed  to  over  850  original  data  points  in  2' -inter¬ 
vals.  Since  the  (14,14)  S.H.  model  corresponds  approximately  to  a  13°-resolution, 
the  h°  data  interval  could  be  replaced  by  a  much  larger  interval  without  affecting 


the  S.H.  adjustment  itself.  For  example,  one  could  use  a  l°-interval  and  fill 
in  the  residuals  at  H°~  or  shorter  intervals  as  needed  for  subsequent  determina¬ 
tions  of  a  detailed  gravity  field.  This  would  be  done  in  the  manner  explained 
below  in  conjunction  with  the  residuals  filled  in  at  all  the  original  data  points 
A  good  resolution  of  the  state  vector  parameters  (six  per  arc)  within  this  S.H. 
adjustment  would  then  require  a  minimum  arc  length  of  6°  allowing  for  7  or  more 
altimeter  observations. 

On  the  other  hand,  the  %°-interval  admits  arcs  as  short  as  3°,  which  is  also 
the  criterion  of  the  minimum  arc  length  from  the  standpoint  of  sufficient  number 
of  intersections  with  other  arcs  (see  the  above  reference).  This  criterion 
allows  for  the  utilization  of  80-90%  of  all  SEASAT  arcs.  During  the  actual  S.H. 
adjustment  of  SEASAT  altimetry  at  AFGL  the  observational  interval  has  been  chosen 
at  H°  and  the  minimum  length  criterion  has  been  adopted  as  3°.  Accordingly,  no 
filling  in  of  residuals  has  been  necessary  for  a  subsequent  P.M.  adjustment  or 
collocation  predictions  resulting  in  a  2°-geoid.  The  S.H.  adjustment  has  en¬ 
compassed  some  6,700  arcs,  of  which  about  one-half  have  been  25°  in  length. 

Since  the  S.H.  adjustment  can  optionally  solve  for  chosen  tidal  parameters 
as  well,  repeating  tracks  covering  different  phase  angles  of  tidal  constituents 
are  particularly  beneficial.  However,  in  the  detection  of  abrupt  sea  surface 
dips  via  the  filled-in  residuals  the  tidal  effects  are  negligible  and  the 
repeating  tracks  lose  much  of  their  usefulness.  Accordingly,  a  good  percentage 
of  arcs  could  be  left  out  of  consideration  when  examining  these  residuals.  It 
is  then  quite  natural  to  eliminate  the  shortest  arcs  whose  state  vectors  are 
the  most  susceptible  to  contamination  by  the  sea  surface  dip  itself.  By  con¬ 
trast,  relatively  long  arcs  crossing  a  trench  area,  such  as  15°  -  25°  arcs,  have 
only  a  small  portion  of  their  length  exposed  to  the  sea  surface  dip;  this  does 
not  allow  for  any  significant  transfer  of  the  negative  anomaly  from  the  residuals 
into  the  state  vector  parameters. 


One  can  thus  accomplish  the  double  task  --  eliminating  both  a  number  of 
unnecessary  arcs  and  most  of  the  arcs  affected  by  sea  surface  dips  themselves  -- 
by  adopting  a  minimum  arc  length  of  15°  (corresponding  to  about  500  original 
data  points)  for  the  purpose  of  this  analysis.  Even  so,  over  two- thirds  of  all 
the  SEASAT  arcs  will  still  be  retained  including  those  along  repeating  tracks. 
This  criterion  could,  of  course,  be  increased  to  20°  or  another  suitable  length. 
Exceptionally,  residuals  on  even  a  relatively  long  arc  can  be  contaminated  in 
the  above  sense  if  its  ground  track  happens  to  coincide  with  the  axis  of  an 
extended  trench.  But  this  is  not  cause  for  concern  by  virtue  of  the  criss¬ 
crossing  pattern  of  SEASAT  tracks.  If,  for  example,  an  ascending  pass  follows 
the  direction  of  a  trench,  a  descending  pass  will  cross  it,  and  will  be  able 
to  detect  the  sea  surface  dip  with  maximum  efficiency. 

The  simplest  approach  to  filling  in  the  residuals  at  data  points  not  used 
in  the  global  S.H.  adjustment  will  now  be  described  in  detail.  The  crucial  role 
in  this  task  is  played  by  the  magnetic  tape  containing  the  original  altimeter 
observations  and  the  initial  state  vector  (s.v.)  parameters  on  each  arc. 

According  to  previous  statements,  during  the  S.H.  adjustment  only  one  out  of  16 
observations  on  this  tape  constitutes  the  input  quantity  leading  to  the  formula¬ 
tion  of  the  pertinent  observation  equation  and,  eventually,  to  the  final  solution 
In  addition  to  the  S.H.  potential  coefficients  and  selected  tidal  parameters,  the 
solution  also  includes  corrections  to  the  s.v.  parameters  on  each  arc.  The 
corrected  s.v.  parameters  can  then  replace  their  initial  values  on  this  tape,  or 
else  should  be  recorded  on  a  separate  tape.  The  key  role  in  the  subsequent  com¬ 
putation  of  the  high-density  residuals  is  played  by  the  original  observations 
coupled  with  the  new  s.v.  parameters.  These  two  kinds  of  quantities,  optionally 
supplemented  by  the  effect  of  tidal  adjustment  (deemed  negligible  in  this 
context),  lead  directly  to  the  values  H„o„  and  N„o„  depicted  in  Fig.  A.  The 
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adjusted  S.H.  coefficients  (stored  on  another  magnetic  tape)  give  rise  to  the 
values  Ha  and  Na.  Thus  the  values  vH  =  -vN  from  (A3.1a-c)  can  be  computed  at 
every  original  data  point. 

In  practice  this  task  can  be  accomplished  with  ease  and  efficiency  upon 
undertaking  what  is  essentially  a  new  iteration  in  the  S.H.  adjustment  process 
in  which,  however,  only  a  few  initial  steps  are  carried  out.  The  highlights  of 
this  process  are: 

a)  all  the  altimeter  observations  are  read  and  utilized; 

b)  the  initial  s.v.  parameters  are  replaced  by  their  adjusted  values; 

c)  the  initial  S.H.  coefficients  are  replaced  by  their  adjusted  values 
(a  similar  step  could  optionally  take  place  with  regard  to  the  tidal 
parameters  as  well); 

d)  the  adjustment  routine  proceeds  to  the  formation  of  the  constant 
terms  of  observation  equations,  stores  them  on  a  magnetic  tape  and 
stops; 

e)  these  terms  are  the  desired  high-density  altimeter  residuals. 

The  last  statement  stems  from  the  fact  that  the  quantities  called  "constant  terms" 
in  conjunction  with  the  initial  values  of  parameters  become  residuals  in  conjunc¬ 
tion  with  the  final  values  of  these  parameters.  As  has  transpired,  the  final 
values  of  the  parameters  have  been  obtained  very  economically  by  utilizing  only 
one-sixteenth  of  the  available  data,  but  serve  subsequently  in  the  evaluation  of 
the  residuals  corresponding  to  the  entire  data  set.  In  spite  of  a  great  number 
of  residuals  to  be  evaluated,  this  step  is  exceedingly  efficient  when  compared 
to  the  complete  adjustment  process. 

Having  now  high-density  and  high-resolution  altimeter  residuals  available, 
one  is  in  a  position  to  detect  abrupt  sea  surface  dips  and  evaluate  their 
significance.  This  can  lead  to  the  detection  of  ocean  trenches  and  other 
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anomalies.  The  altimeter  residuals  should  be  considered  against  the  yardstick 
of  one  siqma  characterizing  the  geoidal  resolution,  such  as  o  presented  in  (A3. 3). 
Clearly,  if  the  magnitude  of  the  (negative)  altimeter  residuals  is  much  larger 
than  over  an  extended  area,  e.g.  over  100  km  of  a  satellite  profile,  a  presence 
of  a  significant  ocean  trench  is  suspected. 

In  practical  terms,  the  residuals  can  be  scanned  for  large  negative  values 
such  as 

VH  <  '2°  • 

If  several  consecutive  residuals  along  a  given  arc  fit  this  characteristic,  the 
pertinent  portion  of  the  arc  can  be  plotted  with  the  rms  residual  also  indicated. 
The  threshhold  number  of  such  consecutive  residuals  can  be  stipulated  as  16, 
corresponding  to  a  i°-segment  of  the  profile.  The  above  two  scrutiny-type  param¬ 
eters,  -2 a  and  i°,  are  subject  to  change  depending  on  the  size  of  trenches  one 
wishes  to  detect.  For  the  detection  of  deep  and  wide  features  of  the  sea  surface 
the  magnitude  of  the  two  parameters  can  be  increased.  Conversely,  for  the  detec¬ 
tion  of  smaller  features  the  magnitude  would  be  lowered.  However,  since  the 
actual  geoidal  variations  themselves  (without  the  contribution  of  trenches)  are 
on  the  order  of  a  with  respect  to  the  adjusted  (14,14)  geoid,  there  is  a  limit 
to  the  size  of  trenches  and  other  anomalies  one  can  detect  with  this  method. 
"Calibration"  tests  using  known  bathymetry,  for  example,  could  furnish  quantita¬ 
tive  information  with  regard  to  this  limitation. 
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